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AN  ASSESSMENT  OF  RESERVOIR  MIXING  PROCESSES 


SECTION  1 .  INTRODUCTION 

Mixing  refers  to  all  physical  transport  processes  and  mechanisms 
which  cause  a  parcel  of  water  and  its  associated  water  quality  constit¬ 
uents  to  blend  with  or  be  diluted  by  another  parcel  of  water.  Examples 
of  specific  transport  processes  that  are  included  in  this  definition  are 
molecular  diffusion,  dispersion,  and  convection,  among  others.  The  rel¬ 
ative  importance  of  any  one  specific  transport  process  (e.g.,  disper¬ 
sion)  will  vary  both  temporally  and  spatially  and  depend  on  the  physical 
characteristics  of  the  water  body,  time  of  year,  and  type  of  forcing 
event  (e.g.,  the  wind). 

In  lakes  and  reservoirs,  mixing  results  from  the  cumulative 
effects  of  external  energy  inputs  such  as  surface  heat  exchange;  absorp¬ 
tion  of  solar  radiation  with  depth;  wind  magnitude  and  direction;  inflow 
magnitude,  density,  and  location;  outflow  magnitude  and  location;  and 
changes  in  project  operation  (i.e.,  withdrawal  depth,  pool  level 
changes,  etc.).  Mixing  is  therefore  dynamic  since  its  effectiveness 
varies  in  response  to  those  dynamic  forcing  events. 

A  thorough  understanding  of  mixing  is  required  since  it  and  the 
resultant  hydrothermal  regime,  which  includes  thermal  stratification,  is 
a  dominant  factor  in  determining  what  takes  place  chemically  and  biolog¬ 
ically  in  a  reservoir.  Sensitivity  analyses  on  reservoir  water  quality 
models  have  shown  the  mixing  coefficient  to  be  highly  sensitive,  indi¬ 
cating  the  importance  of  having  an  accurate  formulation  for  mixing 
(Thornton  et  al.  1979). 

The  objective  of  this  study  is  twofold: 

a.  To  review  and  document  the  major  mixing  mechanisms  in  Corps  of 
Engineers  (CE)  reservoirs  and  their  importance  to  reservoir 
water  quality. 

b.  To  develop  a  mathematical  algorithm  for  one-dimensional  water 
quality  models  (i.e.,  CE-QUAL-R1,  Environmental  Laboratory 
1982)  which  realistically  represents  all  major  mixing  pro¬ 
cesses  occurring  in  CE  reservoirs. 


The  review  will  emphasize  the  transport  and  dispersion  of  sub¬ 
stances  due  to  a  given  hydrothermal  regime  and  not  the  specific  details 
of  the  hydrodynamic  processes  themselves.  The  mathematical  algorithm 
must  be  sufficiently  general  to  be  applicable  to  the  majority  of  CE 
reservoirs  which  vary  in  size  and  location  and  are  therefore  driven  by 
different  hydrometeorologic  forcing  events.  The  algorithm  must  also  be 
able  to  accurately  simulate  changes  in  the  mixing  regime  due  to  changes 
in  project  operations,  such  as  a  change  in  conservation  pool  level  or 
withdrawal  depth. 
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SECTION  2.  BACKGROUND  AND  REVIEW 


2.1  Introduction 

Although  the  ultimate  objective  of  this  review  is  knowledge  of  the 
sum  or  cumulative  effects  of  all  reservoir  mixing  processes  or  mecha- 
nisms  and  their  impact  on  reservoir  water  quality,  it  is  first  necessary 
to  review  density  stratification  and  the  individual  transport  processes. 
In  many  instances,  this  is  not  a  simple  task  because  the  individual 
transport  mechanisms  act  together,  reinforcing  one  another  with  unknown 
synergistic  effects.  In  addition,  mixing  determines  the  observed  ther¬ 
mal  structure,  but  the  thermal  (density)  stratification  modifies  the 
mixing  regime.  It  is  therefore  impossible  to  discuss  reservoir  mixing 
without  knowledge  of  thermal  stratification  and  to  understand  thermal 
stratification  and  reservoir  water  quality  without  an  understanding  of 
the  individual  mixing  processes. 

In  this  section,  the  concepts  of  thermal  stratification  and  poten¬ 
tial  energy  will  first  be  reviewed.  Then,  definitions  for  several 
fundamental  transport  processes  will  be  presented.  This  background 
information  will  be  used  in  the  next  section  to  construct  a  complete 
picture  of  reservoir  mixing  based  on  energy  sources  and  concepts  and  to 
determine  the  significance  of  mixing  on  reservoir  water  quality. 

2.2  Density  Stratification 

2.2.1  Definitions  and  Importance 

Density  stratification  Is  the  nonhomogeneous  layering  of  a  fluid 
due  to  differences  in  density.  In  reservoirs,  density  stratification  is 
primarily  caused  by  temperature  (i.e.,  thermal  stratification),  but  den¬ 
sity  differences  resulting  from  variations  in  suspended  and  dissolved 
solids  concentrations  can  also  be  important.  Density  and/or  thermal 
stratification  therefore  implies  incomplete  vertical  mixing. 

Many  deep  reservoirs  exhibit  the  classical  three-layer  stratifi¬ 
cation  (Figure  la).  The  epilimnion  is  the  upper  stratum  of  warm. 
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turbulent  water.  It  Is  usually  characterized  by  relatively  uniform 
temperatures.  The  deep,  cold,  relatively  undisturbed  region  is  termed 
the  hypolitrnion.  Between  the  epilimnion  and  hypolimnion  is  a  transition 
zone,  the  metalimnion,  which  is  characterized  by  a  strong  temperature 
(density)  gradient.  The  plane  of  inflection  or  of  maximum  temperature 
gradient  is  termed  the  thermocline.  Other  definitions  have  been  pro¬ 
posed  for  the  thermocline  (e.g.,  the  1°  C/m  criterion),  but  they  are  not 
used  herein.  Another  term  illustrated  in  Figure  la  is  the  "mixed 
layer."  The  mixed  layer  is,  as  implied  by  the  name,  the  overlying 
isotropic  layer.  Since  the  mixed  layer  refers  to  the  instantaneous 
depth  of  the  overlying  isotropic  layer,  it  differs  from  the  epilimnion, 
which  is  an  averaged  mixed  layer,  in  two  respects.  First,  the  mixed 
layer  depth  is  usually  less  than  the  depth  of  the  epilimnion.  Second, 
it  is  much  more  dynamic  than  the  epilimnion. 

Many  shallow  reservoirs  with  short  hydraulic  residence  times  do 
not  exhibit  a  classical  three-layered  system  (e.g..  Figure  lb).  Instead 
of  having  a  well-defined  hypolimnion,  the  metalimnion  appears  to  extend 
to  the  reservoir  bottom.  It  is  also  possible  to  have  systems  that  do 
not  appear  to  have  a  well-defined  epilimnion. 

In  cold  regions,  reservoirs  may  inversely  stratify  during  the 
winter  months  (Figure  lc) .  At  low  temperatures  (near  4°  C) ,  even  small 
quantities  of  solids  can  significantly  alter  the  water  density  and 
impact  the  observed  thermal  structure.  In  general,  there  is  no  well- 
defined  epilimnion  and  metalimnion  in  inversely  stratified  systems.  The 
thickness  of  the  zone  of  density  gradient  (i.e.,  the  epilimnion  and 
metalimnion)  under  the  ice  varies  from  a  metre  or  two  in  small  reser¬ 
voirs  to  tens  of  metres  in  large,  deep  reservoirs.  In  deep  reservoirs, 
the  effects  of  hydrostatic  pressure  (i.e.,  depth)  also  modify  the  den¬ 
sity  distribution  (Farmer  and  Carmack  1982) . 

Many  other  terms  are  used  to  describe  the  thermal  structure  of  a 
lake.  These  are  defined  in  Hutchinson  (1957),  Ruttner  (1963),  Wetzel 
(1983),  among  others.  In  general,  these  terms  are  related  either  to  the 
number  of  turnovers  (i.e.,  periods  of  complete  vertical  mixing)  occur¬ 
ring  within  a  lake  or  to  the  strength  of  the  stratification.  It  is. 
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however,  important  to  distinguish  between  holomictic  and  meromictic 
lakes.  In  holomictic  lakes,  the  entire  water  column  completely  circu¬ 
lates  or  turns  over.  Lakes  that  cannot  circulate  completely  and  exhibit 
a  deep  stratum  that  is  perennially  stagnant  in  the  water  column  are 
termed  meromictic  lakes  (Hutchinson  1957).  The  reason  for  this  condi¬ 
tion  may  be  physical,  chemical,  or  biological.  Turbidity- Induced 
meromixls  has  been  observed  in  a  CE  reservoir  (Larson  1979).  Whatever 
the  cause  of  meromixis,  it  can  have  a  profound  effect  on  the  temperature 
and  mixing  structure  and,  consequentially,  on  the  water  quality  of  lake. 
Bottom  withdrawal  car  be  an  effective  management  alternative  to  elimi¬ 
nate  this  undesirable  condition. 

Thermal  stratification  is  important  because  all  chemical  and  bio¬ 
logical  processes  are,  to  some  extent,  temperature  dependent.  More 
important,  however,  is  the  layering  of  the  lake  which  can  isolate  the 
metalimnion  and  hypolimnion  from  light  and  transfers  across  the  air / 
water  interface,  resulting  in  vertical  variations  in  water  quality. 

2,2.2  Factors  Affecting  Stratification 

The  principal  factors  Influencing  the  formation,  strength,  and 
extent  of  stratification  are  the  density  of  water;  solar  radiation  and 
the  heat  transfer  at  the  air /water  interface;  and  the  mixing  resulting 
fri.’’,  advection  (inflows  and  outflows)  and  wind-induced  phenomena. 

It  is  well  known  that  the  density  of  water  varies  with  temperature 
(Figure  2) ,  The  importance  of  this  variation  in  determining  the  distri¬ 
bution  of  heat  within  a  lake  was  originally  documented  by  Birge  (1910). 
Two  factors  are  important.  First,  the  maximum  density  of  water  occurs 
at  approximately  4°  C.  Second,  water  density  decreases  at  an  escalating 
rate  with  both  increasing  and  decreasing  temperatures  from  4°  C.  There 
is  over  an  order  of  magnitude  difference  in  the  energy  requirements  to 
mix  or  destratify  a  1°  C  temperature  difference  at  25°  C  than  at  5°  C. 

The  energy  available  to  warm  the  waters  of  a  reservoir  ultimately 
comes  from  solar  radiation,  which  varies  seasonally  and  with  latitude 
(Figure  3).  The  seasonal  variation  of  solar  radiation  follows  a  sinu¬ 
soidal  curve  with  a  maximum  in  late  June.  In  addition,  diurnal  cycles 
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Figure  2.  Variations  of  water 
density  with  temperature 


Figure  3.  Seasonal  variations  in 
solar  radiation  (after  Koberg  1962) 

also  occur.  Water  temperatures  respond  to  both  of  these  cycles  with  a 
slight  delay.  Solar  radiation  is  absorbed  at  the  water  surface  and 
selectively  with  depth  depending  on  the  wavelength  of  the  light,  prop¬ 
erties  of  the  water,  and  the  matter  suspended  in  the  water.  This 
absorption  is  usually  assumed  to  be  exponential  with  depth  (i.e.,  Beer's 
Law),  but  surface  effects  result  in  minor  discrepancies  in  the  top  metre 
or  so  of  a  lake  (Figure  4) . 

In  contrast  to  heating,  cooling  of  a  water  body  can  occur  only  at 
the  water  surface.  It  is  therefore  possible  for  the  surface  water  to 
decrease  in  temperature  while  deeper  water  increases  in  temperature.  If 
the  temperature  of  the  surface  water  drops  below  the  temperature  of  the 
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deeper  water  and  still  remains  above  C,  the  water  column  becomes 
thermally  unstable  and  natural  convection  and  mixing  commence. 

Mixing  causes  the  shape  of  the  vertical  temperature  profile  to 
change  from  an  exponential  decrease  with  depth  (i.e.,  warming  due  to 
absorption  of  solar  radiation)  to  the  classical  profile  of  a  well-mixed 
layer  overlying  a  zone  of  temperature  gradient  (Figure  5) .  Reservoir 
mixing  is  discussed  in  detail  in  Section  3. 

2,2,3  Temporal  Variations 

Because  the  factors  which  determine  the  thermal  stratification  are 
always  changing,  the  thermal  structure  is  always  undergoing  change. 

There  are,  however,  three  distinct  cycles  of  importance;  annual,  synop¬ 
tic,  and  diel.  The  annual  cycle  has  a  period  of  365  days  and  exhibits 
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b.  Solar  heating 
plus  wind 


Figure  5.  Effect  of  mixing  on  the  thermal  profile 


seasonal  changes  in  temperature  resulting  from  seasonal  changes  in  solar 
radiation,  air  temperature,  wind,  and  flow.  Synoptic  cycles  typically 
have  periods  of  5  to  7  days  and  correspond  with  the  passage  of  major 
weather  systems  (i.e.,  warm  or  cold  fronts).  Diel  cycles  have  a  period 
of  2 4  hr  and  correspond  to  daytime  heating  and  nighttime  cooling. 

Annual  cycle.  The  annual  temperature  cycle  for  a  large  reservoir 
with  long  residence  time  is  shown  in  Figure  6,  During  the  period  of 
spring  turnover  (late  February  to  early  March),  the  entire  water  column 
undergoes  complete  mixing.  The  density  differences  at  the  low  tempera¬ 
tures  (4°  to  5®  C)  are  insufficient  to  prevent  complete  mixing.  As  the 
water  column  warmed,  density  differences  increase  and  it  becomes  more 
difficult  to  mix  the  entire  water  column.  Stratification  started  to 
form  near  the  bottom  of  the  lake  (late  March)  because  the  density  dif¬ 
ferences  and  resulting  buoyancy  forces  were  small  compared  to  the 
kinetic  energy  input  (i.e.,  wind).  As  the  solar  radiation  increased, 
water  temperatures  increased,  density  difference  increased,  and  the 
thermocline  moved  upward  because  the  kinetic  energy  input  could  not 
overcome  the  ever-increasing  buoyancy  forces  of  density  stratification. 
The  minimum  thermocline  depth  was  achieved  about  the  time  of  summer 
solstice  or  time  of  maximum  heat  input  (late  June).  Once  stratification 
formed,  the  hypolimnetic  temperatures  remained  relatively  constant  until 
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Figure  6.  Annual  temperature  variations  in  a 
deep  reservoir  with  a  long  residence  time 

the  surface  temperatures  cooled  to  near  the  hypolimnetic  temperatures 
and  complete  vertical  mixing  of  the  water  column  occurred  (i.e.,  fall 
overturn) . 

The  increased  mixing  resulting  from  inflows  and  outflows  in  reser¬ 
voirs  can  result  in  deviations  from  the  above  example.  For  instance, 
periods  of  overturn  may  be  extended  and  become  more  frequent,  and  the 
slope  of  the  isotherms  may  be  increased.  If  bottom  withdrawal  is  used, 
hypolimnetic  temperatures  may  increase.  For  example,  in  Figure  7,  the 
hypolimnetic  temperatures  in  Beltzville  Lake,  Pennsylvania, during  1976 
remained  relatively  constant  when  only  small  amounts  of  water  were 
released  through  the  lower  ports.  In  contrast,  during  1972,  when  large 
quantities  of  water  were  released  through  the  floodgates  in  response  to 
Hurricane  Agnes,  the  hypolimnetic  temperatures  increased.  Other 
variations  such  as  pumped-storage  operations  can  also  increase  mixing 
and  hypolimnion  temperatures  (Figure  8) . 


BELTZVILLE  LAKE 


Synoptic  variations.  Examples  of  synoptic  temperature  variations 
are  shown  in  Figure  9.  Temporary  periods  of  stratification  occur  during 
periods  of  warm,  calm  weather  and  are  destroyed  during  periods  of  cold, 
windy  weather.  Synoptic  variations  are  on  the  order  of  a  few  degrees 
Celsius. 

Dlel  variations.  Diel  variations  are  typically  on  the  order  of  1° 
to  2°  C  but  can  be  as  large  as  7°  C  or  more.  The  actual  magnitude  will 
depend  on  the  depth  of  the  upper  mixed  layer,  the  amount  of  surface  mix¬ 
ing,  and  the  quantity  of  solar  insolation.  The  deeper  the  mixed  layer 
and/or  the  larger  the  surface  mixing,  the  smaller  the  diurnal  variation. 
A  typical  diel  variation  is  shown  in  Figure  10. 

2.2.4  Horizontal  Variations 

Horizontal  variations  in  temperature  and  stratification  occur  as  a 
result  of  differential  heating,  inflow,  or  mixing.  Differential  heating 
occurs  when  the  smaller  volume  of  water  in  the  coves,  littoral  zones, 
and  headwaters  of  an  impoundment  warms  or  cools  more  rapidly  than  the 
open-water  regions.  In  large  lakes,  this  phenomenon  is  significant  and 
results  in  the  formation  of  thermal  bars.  Similarly,  rivers  flowing 
into  a  reservoir  may  be  of  different  temperatures,  creating  longitudinal 
variations.  Horizontal  variations  resulting  from  a  river  inflow  are 
illustrated  in  Figure  11.  Horizontal  variations  resulting  from  a  river 
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Figure  8.  Effect  of  pumped-storage  operations  on  the  thermal 
structure  of  Kinzua  Lake,  Pennsylvania  (after  Dortch  1981) 

(to  convert  feet  to  metres,  multiply  by  0.3048) 


inflow  typically  create  temperature  differences  of  1°  or  2°  C  or  more. 
Spatial  variations  in  both  the  horizontal  and  vertical  directions  can 
also  occur  as  a  result  of  seiching  and  upwelling  (see  Sections  3,2.2  and 
3.3.2).  These  variations  are  highly  dynamic. 
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Figure  9.  Example  of  synoptic  variations  in  water  temperature 
DeGray  Lake,  Arkansas 


ELEVATION.  M 


TEMPERATURE,  °C 


JULY 

Figure  10,  Example  of  a  typical  diel 
temperature  cycle,  DeGray  Lake,  Arkansas 


Figure  11.  Horizontal  variations  in  temperature  resulting 
from  river  inflow,  DeGray  Lake,  Arkansas,  2  Apr  74 
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2.2.5  Data  Interpretation  and  Generalizations 

Temperature  data  interpretation.  Much  of  our  knowledge  concerning 
mixing  in  reservoirs  and  lakes  is  based  on  the  interpretation  of  verti¬ 
cal  temperature  profile  data  (i.e.,  stratification  data).  In  addition 
to  understanding  the  factors  which  affect  stratification  and  its  tempo¬ 
ral  and  spatial  variations  (Sections  2. 2. 2-2. 2. 4) ,  the  correct  interpre¬ 
tation  of  stratification  data  llso  requires  detailed  knowledge  of  the 
field  sampling  (e.g.,  station  location,  sampling  time,  etc.),  hydro¬ 
meteorological  conditions  (e.g.,  inflows,  wind  speed  and  direction), 
project  operation  (e.g.,  release  rates  and  outlet  locations),  and  reser¬ 
voir  mixing  processes.  There  should  be  a  logical,  physically  based 
explanation  for  all  observations  (Ford  1978) . 

For  example,  most  field  data  are  collected  during  daylight  hours 
and  may  represent  maximum  daily  water  temperature  and  minimum  mixed 
layer  depths  (see  Figure  10).  If  there  is  wind,  warm  surface  water  may 
be  pushed  to  one  side  of  the  reservoir  (Figure  12)  and/or  the  thermo- 
cline  may  tilt,  changing  the  depth  of  the  upper  mixed  layer. 

When  interpreting  stratification  data,  it  is  usually  beneficial  to 
draw  isotherms  relative  to  the  water  surface  (e.g.,  Figure  13).  For 
reservoirs,  it  is  essential  that  the  variation  in  water  surface  eleva¬ 
tion  be  considered  since  there  is  the  potential  for  large  variations  in 
water  levels  and  since  the  slope  of  the  isotherms  indicates  the  degree 
of  mixing.  The  relatively  flat  isotherms  in  the  lower  metalimnion 
(i.e.,  16°  and  18°  C)  during  midsummer  indicate  little  mixing  in  this 
region  while  the  steep  slope  of  the  22°  and  24°  C  isotherms  in  the  upper 
metalimnion  during  the  same  period  possibly  indicates  more  intensive 
mixing.  If  the  variation  in  the  water  surface  elevation  is  not  taken 
into  consideration,  the  slope  of  the  isotherms  may  be  misrepresented. 
Short-term  variations  In  the  isotherms  (i.e.,  order  of  days)  indicate 
seiching  and/or  Internal  waves  (Section  3.2,2)  and  should  be  averaged 
when  comparing  isotherm  slopes.  This  should  not  be  a  problem  if  the 
time  interval  between  sampling  dates  is  greater  than  2  weeks. 

Once  the  isotherms  are  constructed,  the  importance  of  inflows, 
outflows,  and  light  penetration  can  be  determined  by  comparing  periods 


WIND 


Figure  12.  Effect  of  wind  on  temperature 
profiles,  McCarrons  Lake,  MiTnesota,  1974 

and  regions  of  mixing  with  the  inflow  quantity  and  placement  (based  on 
density  or  temperature)  (Figure  14a) ,  the  outflow  quantity  and  depth 
(Figure  14b) ,  and  the  seasonal  variation  in  twice  the  Secchi  disc  depth 
(i.e.,  ca.  1  percent  light  level)  (Figure  14c).  For  example,  the  inten¬ 
sive  mixing  that  occurred  during  April  and  May,  as  evidenced  by  the 
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Figure  13.  Construction  of  isotherms  relative  to 
water  surface  level,  DeGray  Lake,  Arkansas,  1978 

steep  slope  of  the  10°  to  18°  C  isotherms,  is  probably  the  result  of  the 
large  outflows  and  change  in  withdrawal  depth  (Figure  14b)  since  light 
(thermal  energy)  did  not  penetrate  to  these  depths  (Figure  14c) .  In 
contrast,  in  Figure  15  the  steep  gradient  of  the  14°  to  16°  C  isotherms 
during  early  June  is  the  result  of  internal  absorption  of  solar  radia¬ 
tion  since  the  light  penetration  increased  significantly  during  this 
period  and  winds  and  flows  were  low. 

One  aspect  of  data  interpretation  that  is  commonly  overlooked  is 
the  error  associated  with  field  measurements.  Because  of  the  intrinsic 
variability  of  temperature  in  reservoirs  and  calibration  error  associ¬ 
ated  with  field  equipment,  historical  temperature  data  should  be 
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Figure  14.  DeGray  Lake,  Arkansas,  1979 
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Figure  15.  Comparison  of  isotherms 
and  photic  zone.  Lake  Calhoun, 
Minnesota,  1974 


considered  no  better  than  ±1°  C.  For  example,  pockets  of  cold  water  in 
the  hypolimnion  of  a  lake  (Figure  16)  are  sometimes  attributed  to  an 
influx  of  cold  ground  water  rather  than  a  data  error.  Such  a  situation 
is  highly  unlikely  since  ground-water  movement  is  slow.  If  it  should 
occur,  it  would  be  apparent  in  the  daily  water  budget  (i.e.,  change  in 
pool  level  or  outflow  rate)  as  well  as  in  the  temperature  profile. 

Generalizations  concerning  stratification.  Several  generaliza¬ 
tions  concerning  stratification  in  CE  reservoirs  can  be  made  based  upon 
a  review  of  temperature  data  collected  at  reservoirs  throughout  the 
United  States.  First  and  foremost,  all  reservoirs  stratify,  albeit  some 
for  only  a  few  hours  and  a  few  degrees  Celsius.  Second,  if  the  mean 
annual  theoretical  hydraulic  residence  time  (i.e.,  volume/mean  annual 
flow  for  period  of  record)  is  greater  than  6  months,  the  stratification 
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Figure  16.  Isotherms  for  DeGray  Lake, 
Arkansas,  1980 


is  dominated  by  meteorological  forcing.  Stratification  in  this  type  of 
lake  is  characterized  by  several  features: 

a.  Since  the  major  factors  responsible  for  stratification  (e.g., 
wind  and  solar  radiation)  act  at  the  air/water  interface, 
horizontal  variations  are  minimized. 

b.  There  is  not  much  variability  in  stratification  from  year  to 
year.  Stratification  forms  and  fall  overturn  occurs  at 
approximately  the  same  time  each  year.  Hypolimnetie 
temperatures  and  thermocline  depths  are  similar  from  year  to 
year  (e.g.,  Table  1). 

£.  The  larger  the  surface  area,  the  more  wind  (kinetic)  energy 
input,  the  longer  the  periods  of  turnover,  the  deeper  the 
upper  mixed  layer. 

d.  The  deeper  the  lake,  the  less  wind  (kinetic)  energy  per  unit 
volume  available  for  mixing,  and  the  stronger  the  stratifi¬ 
cation. 

e.  Lakes  in  similar  geographic  areas  will  be  exposed  to  similar 
hydrometeorological  conditions  and  exhibit  similar  stratifi¬ 
cation  patterns.  For  example.  Lakes  DeGray,  Greeson,  and 
Ouachita  are  located  within  55  km  of  each  other  in  the 
Ouachita  Mountains  of  southwestern  Arkansas,  are  deep  (i.e., 
>30  m) ,  have  residence  times  greater  than  12  months,  and 
therefore  exhibit  similar  stratification  profiles 

(Figure  17). 

Third,  if  the  mean  annual  theoretical  hydraulic  residence  time  is 
small  (less  than  20  days) ,  the  system  is  advectively  (inflows  and 
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Table  1 

Variation  in  August  Stratification  Characteristics 
DeGray  Lake,  Arkansas 


Year 

Withdrawal 

Location 

Hypolimnetic 
Temperature,  °C 

Thermo cline 
Depth,  m 

1974 

Surface 

7.0 

8 

1975 

Surface 

8.0 

8 

1976 

Surface 

7.0 

7 

1977 

Surface 

6.0 

7 

1978 

Surface 

6.0 

7 

1979 

Lower 

6.5 

8 

1980 

Lower 

7.5 

8 

1981 

Lower 

6.8 

10 

1982 

Lower 

6.2 

8 

1983 

Surface 

7.0 

8 

outflows)  dominated,  is  weakly  and  intermittently  stratified,  and  will 
exhibit  characteristics  similar  to  the  inflows.  During  storms,  the 
residence  times  of  these  projects  may  be  only  a  few  hours,  indicating 
that  detailed  knowledge  of  the  storm  runoff  may  be  required  to  explain 
and  predict  the  density  (thermal)  structure  of  these  projects.  Fourth, 
if  the  mean  annual  theoretical  hydraulic  residence  time  is  greater  than 
20  days  but  less  than  6  months,  the  stratification  will  be  controlled 
both  by  meteorological  forcing  and  advection.  In  these  projects,  it  is 
possible  for  the  stratification  pattern  to  vary  significantly  from  year 
to  year. 

Specific  guidance  concerning  criteria  to  evaluate  the  stratifica¬ 
tion  potential  of  a  reservoir  is  given  in  Appendix  A. 

2.3  Potential  Energy 

Closely  related  to  density  (thermal)  stratification  is  the  concept 
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Figure  17.  Comparison  of  temperature 
profiles  from  Lakes  DeGray,  Greeson, 
and  Quachita,  Arkansas,  1962 


Temperature,  *C 


of  potential  energy  (PE*) .  Potential  energy  is  a  form  of  stored  energy 
that  a  system  possesses  because  of  its  configuration.  The  PE  associated 
with  density  stratification  is  dependent  on  gravity  and  conservative 
force,  and  can  be  fully  recovered  and  converted  into  kinetic  energy. 

The  FE  of  a  reservoir  is  defined  as: 


(1) 


PE 


mgH 


8  z  A(z) 


p(z,t)  dz 


where 

m  ■  the  total  mass  of  the  reservoir,  kg 

2 

g  -  the  acceleration  due  to  gravity,  m/sec 

H  *  height  of  the  center  of  mass  of  the  reservoir,  m 

Z  *  the  maximum  elevation,  m 


*  For  convenience,  symbols  and  abbreviations  are  listed  in  the  Notation 
(Appendix  B) . 
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z  *  the  elevation  above  the  reservoir  bottom,  m 

A(z)  -  the  horizontal  area  of  the  reservoir  at  elevation  z  , 

2 

m 

p(z,t)  ■  the  reservoir  density  at  elevation  z  and  time  t  , 
kg/m3 

The  PE  of  a  reservoir  can  therefore  be  changed  by  heating  or  cooling  to 
modify  the  water  density  and/or  by  changing  the  elevation  of  the  center 
of  mass.  For  a  well-mixed  body  of  water,  the  center  of  mass  is  the  cen¬ 
ter  of  volume.  For  a  stratified  body  of  water,  the  surface  waters  are 
less  dense  and  the  center  of  mass  is  deeper.  The  potential  energy  con¬ 
cept  is  similar  to  Birge's  wind  work  and  Schmidt's  stability 
(Hutchinson  1957).  For  example,  Schmidt's  stability  is  the  change  in  PE 
(using  Equation  1)  between  an  initial  stratified  density  profile  and  the 
resulting  isothermal  condition  following  complete  mixing. 

The  significance  of  the  potential  energy  concept  is  its  relation¬ 
ship  with  work  and  kinetic  energy  (KE)  as  defined  in  classical  physics. 
It  indicates  that  work  is  required  to  mix  a  stably  stratified  fluid 
since  the  center  of  mass  must  be  raised  against  the  force  of  gravity. 

For  example,  completely  mixing  the  two-layered  fluid  in  Figure  18 
changed  the  PE  by  ApVH/8  and  raised  the  center  of  mass  to  H/2  . 

The  PE,  as  defined  in  Equation  1,  has  limitations  when  applied  to 
reservoirs  because  the  horizontal  areas,  A(z) ,  are  much  larger  and  domi¬ 
nate  over  the  small  differences  in  water  density,  p(z,t).  Additionally, 
the  PE  as  defined  by  Equation  1  decreases  as  stratification  Increases. 

A  more  practical  formulation  that  increases  as  stratification  increases 
is  the  relative  potential  energy  (RPE) : 

Z 

rm 

RPE  “J  g  z  A(z)  pm  -  p(z,t)  dz  (2) 

where  p  ■  maximum  density  in  water  column  . 

If  p  is  a  constant,  the  magnitude  of  the  change  in  PE  between  two 
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Figure  18.  Change  in  potential  energy 
resulting  from  destratification  of  a 
two-layered  system 

states  computed  from  either  Equation  1  or  Equation  2  is  the  same  (Ford 
1976)  but  the  signs  are  opposite. 

The  RPE's  for  DeGray  and  Ouachita  Lakes,  Arkansas,  are  shown  in 
Figure  19.  In  February,  the  reservoirs  are  isothermal  and  the  entire 
water  column  is  mixed  vertically.  The  center  of  mass  is  therefore 
located  at  the  center  of  volume,  and  the  RPE  is  small  since  the  water 
density  of  the  entire  column  is  near  .  The  maximum  RPE  occurs  in 
mid-  to  late-July  and  coincides  with  the  time  of  maximum  stratification. 
As  the  reservoirs  cool,  mixing  occurs,  and  the  RPE  decreases. 

Figure  19  also  illustrates  two  features  of  Equation  2.  First,  the 
larger  the  lake,  the  larger  the  horizontal  areas,  A(z),  and  the  larger 
the  RPE.  In  Figure  19,  Lake  Ouachita  is  significantly  larger  than 
DeGray  Lake.  Second,  the  weaker  the  stratification,  the  smaller  the 
RPE.  In  1976,  DeGray  Lake  was  operated  with  surface  withdrawal  and  had 
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a  stronger  stratification  than  in  1979  when  it  was  operated  with  bottom 
withdrawal.  The  RPE  was  therefore  greater  for  1976  than  for  1979. 

To  illustrate  how  the  RPE  can  be  used  to  explain  the  annual 
thermal  stratification  cycle,  the  RPE  for  DeGray  Lake  in  1979  is 
compared  with  the  turbulent  kinetic  energy  (TKE)  input  from  the  wind  and 
inflows  in  Figure  20.  Specifics  concerning  the  computation  of  TKE  from 
the  wind  and  inflow  rates  are  described  in  Sections  3.2.2  and  3.2.5, 
respectively. 

In  February,  the  reservoir  was  nearly  isothermal.  The  RPE  waG 
small  since  water  densities  were  nearly  constant.  With  little  or  no 
stratification,  the  RPE  was  not  large  enough  to  prevent  complete  mixing 
by  the  IKE.  As  the  water  column  warmed,  stratification  increased,  RPE 
increased,  and  it  became  more  difficult  for  the  TKE  to  mix  the  entire 
water  column.  Stratification  started  to  form  in  March  at  the  bottom  of 
the  lake  (Figure,  6)  because  the  density  differences  and  the  resulting 
RPE  were  small  compared  to  the  TKE  input  from  the  wind  (Ford  and  Stefan 
1980a).  As  the  solar  radiation  increased,  water  temperatures  increased 
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and  the  density  differences  increased;  the  thermocline  moved  upward 
because  the  TKE  input  could  not  overcome  the  ever-increasing  RPE  (i.e., 
buoyancy  forces).  The  maximum  RPE  occurred  in  mid-  to  late-July  and 
coincided  with  the  time  of  maximum  stratification.  Once  stratification 
formed,  the  hypolimnion  temperature  remained  relatively  constant  until 
fall  overturn.  As  the  reservoir  cooled,  stratification  decreased  and 
the  RPE  decreased,  allowing  the  available  TKE  to  deepen  the  thermocline 
further,  and  eventually  resulted  in  complete  vertical  mixing.  When 
comparing  the  RPE  of  stratification  with  the  TKE  input  it  is  important 
to  note  that  the  RPE  varies  significantly  over  the  seasonal  time  scale 
(i.e.,  the  time  scale  of  stratification)  while  the  TKE  is  introduced  at 
much  shorter  time  scales  (i.e.,  passage  of  storms). 

Figure  21  compares  the  TKE  from  the  wind  with  the  TKE  from  the 
inflows.  Except  in  May  (the  wettest  month),  the  TKE  from  wind  is 
greater,  indicating  its  dominance  over  advection  (inflows)  as  a  source 
of  energy  and  as  a  mechanism  for  mixing. 


29 


2.0 


l.S 


®  1.0 


0 


— 

— 

1 - 

1 - 

1 - 

T 

— 

l 

1 

+ 

♦ 

TKE 

FROM 

WIND 

+ 

© 

TKE 

FROM 

INFLOWS 

+ 

• 

+ 

+ 

+  ' 

- 

o 

<9 

© 

■f 

+ 

* 

+ 

- 

o 

■ 

1 

1 

_L 

© 

1  A 

I  © 

A 

* 

© 

L 

JAN 

FEB 

MAR 

APR 

MAY 

JUN 

JUL 

AUG 

SEP 

OCT 

NOV 

DEC 

Figure  21.  Comparison  of  TKE  from  the  wind  and  inflows, 
DeGray  Lake,  Arkansas,  1979 


2.4  General  Transport  Processes 

There  are  several  fundamental  transport  or  mixing  processes  which 
cause  a  parcel  of  fluid  to  blend  with  other  parcels  and  which  occur  in 
any  fluid  system,  including  reservoirs.  These  are  advection,  shear, 
molecular  diffusion,  turbulence  and  turbulent  diffusion,  entrainment, 
convection,  and  dispersion.  Each  of  these  processes  will  be  defined  to 
form  a  common  base  for  the  remainder  of  the  report.  The  nomenclature  of 
Fischer  et  al.  (1979)  will  be  followed. 


2.4.1  Advection 

Advection  is  transport  by  the  mean  motion  of  the  fluid 

(Figure  22a).  If  the  fluid  is  moving  with  velocity  U  ,  which  has 

components  of  U  ,  U  ,  and  U  in  the  x,  y,  z  directions, 
x  y  z 

respectively,  then  the  advective  transport  in  the  direction  of  the 
velocity  vector  U  is  UC  where  C  is  the  concentration  of  the 
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solute.  The  advective  transport  in  the  x,  y,  and  z  directions  is, 

therefore,  U  C  ,  U  C  ,  and  U  C  ,  respectively.  Since  the  units  of  U 

x  y  z  3  3 

are  L/t  (length/tioe)  and  the  units  of  C  are  M/L  (mass/length  ), 

2 

the  units  of  UC  are  M/L  / t,  which  represents  the  rate  of  mass 

2 

transport  through  the  unit  area  L  .  For  transport  by  advection  only, 
there  will  be  no  change  in  the  concentration  C  provided  no  other 
transport  mechanisms  act  on  the  parcel  of  fluid. 

To  be  transported  by  advection,  a  parcel  of  fluid  must  be  acted 
upon  by  a  force.  For  rivers,  this  force  is  gravity.  In  reservoirs, 
advection  may  be  caused  by  inflows,  outflows,  and  wind  shear  at  the 
air/water  interface. 

2.4.2  Shear 

Shear  is  advection  at  different  speeds  at  different  locations 
(i.e.,  flows  with  velocity  gradients  are  shear  flows)  (Figure  22c). 
Friction  or  shear  forces  at  boundaries  cause  velocities  to  be  less  near 
the  boundaries  than  in  the  center  of  the  flow  field.  Shear  flows  vary 
from  simple  logarithmic  profiles  typically  found  in  open-channel  flow  to 
complex  flow  patterns  in  lakes  where  velocity  vectors  vary  temporally 
and  spatially  both  in  magnitude  and  direction.  In  shear  flows,  mixing 
or  spreading  in  the  direction  of  flow  Is  caused  primarily  by  velocity 
gradients.  Since  the  variability  in  the  velocity  profile  is  propor¬ 
tional  to  the  shear  velocity  w*  (w^  >/r/p ,  where  x  *=  shear  stress  and 
p  =  density)  and  independent  of  the  mean  velocity  U  ,  the  longitudinal 
mixing  in  a  shear  flow  depends  on  the  shear  velocity  not  the  mean  flow 
velocity  (Fischer  et  al.  1979). 

2.4.3  Molecular  Diffusion 

Molecular  diffusion  is  a  process  by  which  a  certain  property  of  a 
fluid  is  transferred  down  a  concentration  gradient  by  the  random  motion 
of  molecules  without  any  overall  transport  of  the  fluid  taking  place. 

It  is  important  to  understand  molecular  diffusion  because  molecular 
diffusion  theory  is  the  basis  for  the  turbulent  diffusion  theory  dis¬ 
cussed  in  Section  2.4.4. 
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The  flux  or  transport  of  mass  or  heat  per  unit  cross-sectional 
area  per  unit  time  in  a  given  direction  by  molecular  diffusion  is; 


q 


m 


-K 

m  9x 


(3) 


where 

7 

K  *  molecular  diffusivity  or  diffusion  coefficient,  m'/sec 
3C 

—  -  concentration  gradient,  mg/Jt/m 

cX 

The  minus  sign  is  required  because  the  transport  is  from  high  to  low 
concentrations.  The  molecular  diffusivity  is  a  function  both  of  the 
solute  and  the  solvent  in  which  it  is  dissolved  (Table  2). 


Table  2 

Molecular  Diffusion  Coefficient  of  Solutes  in  Water  at  20°  C 


Substance 

Diffusion  Coefficient 

Temperature 

1.42  x  iff7 

Dissolved  oxygen 

1.80  x  10-9 

Carbon  dioxide 

1.77  x  10"9 

Nitrogen 

1.64  x  10~9 

Sodium  chloride 

1.35  x  10*9 

SOURCE:  "CRC  Handbook  of  Tables  for  Applied  Engineering  Science" 
(Bolz  and  Tuve  1976). 

The  diffusion  equation: 


9C  _  9^C 

9t  m  „  2 
9x 


w 


describes  how  mass  is  transferred  by  molecular  or  Fickian  diffusion. 
The  solution  of  this  equation  for  an  initial  slug  of  mass  M  released 
at  a  point  results  in  a  Gaussian  distribution  of  concentration  versus 
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distance  at  a  fixed  time.  Two  properties  of  this  solution  r^a  impor¬ 
tant.  First,  differences  in  mean  concentrations  are  always  reduced. 
Second,  the  variance  always  grows  linearly  with  time  (Figure  22b). 

2.4.4  Turbulence  and  Turbulent  Diffusion 

Turbulence  is  the  most  important  source  for  mixing  in  lakes 
(Ottesen  Hansen  1978).  It  is  therefore  imperative  that  turbulence  be 
understood  before  mixing  in  lakes  (or  reservoirs)  can  be  understood. 

Most  flows  occurring  in  nature  (e.g.,  atmospheric  and  surface 
water  flows)  are  turbulent.  Two  notable  characteristics  of  these  flows 
are:  (a)  velocities  and  concentrations  at  a  point  in  a  turbulent  flow 

are  unsteady,  and  (b)  mixing  is  much  faster  in  turbulent  flow  than  in 
laminar  flow. 

Turbulence  is  sometimes  described  as  a  family  of  eddies  (i.e., 
rotating  regions  of  fluid).  These  eddies  or  scales  of  turbulence  can 
range  in  size  from  the  physical  limits  of  the  flow  (i.e.»  physical 
dimensions  of  a  reservoir)  down  to  molecular  motion.  At  the  smallest 
scales  of  turbulence,  viscosity  is  important  and  the  motion  dissipates 
into  heat.  One  disadvantage  of  portraying  turbulence  as  a  family  of 
eddies  is  that  it  is  difficult  to  separate  wave  motions  from  turbulence. 
According  to  Steward  (1959),  a  more  precise  definition  of  turbulence  is: 

A  fluid  is  said  to  be  turbulent  if  each  component  of  the 
vorticity  is  distributed  irregularly  and  aperiodically  in 
time  and  space,  if  the  flow  is  characterized  by  a  transfer 
of  energy  from  larger  to  smaller  scales  of  motion,  and  if 
the  mean  separation  of  neighboring  fluid  particles  tends  to 
increase  with  time. 

Turbulent  flows  are  therefore  irregular  (random),  diffusive  (pro¬ 
duce  mixing) ,  rotational  (three-dimensional  vorticity  fluctuations) , 
time  varying,  and  dissipative  (decay  rapidly  without  a  continual  source 
of  energy)  Turbulent  flows  are  characterized  by  large  Reynolds  numbers 
(order  of  10^  in  lakes)  where  energy  propagates  slowly  with  the  speed  of 
the  fluid  motion.  In  contrast,  waves  can  distort  a  density  distribution 


34 


but  cannot  permanently  change  the  stratification  profile  unless  the 
waves  break  to  produce  mixing;  waves  are  dispersive  but  not  dissipative, 
and  energy  is  transferred  rapidly  through  the  fluid  (Turner  1973) . 

The  energy  associated  with  turbulence  or  TKE  is  defined  as 

TKE  -  (u2  +  u2  +  u2)/2  (5) 

x  y  z 

where  u  ,  u  ,  u  *=  turbulent  fluctuations  in  the  x,  y,  and  z  velocity 
x  y  z 

components,  m/sec. 

The  fluctuating  velocity  components  are  related  to  the  instantaneous  and 
mean  velocity  components  by 


u  -  u 

+  u 

(6a) 

x  X 

X 

u  -  u 

+  u 

(6b) 

y  y 

y 

u  -  u 

z  z 

+  u 

z 

(6c) 

U  ,  U  ,  U  -  instantaneous  velocity  components,  m/sec 
x  y  z 

Ux,  U  ,  Uz  =  mean  velocity  components,  m/sec 

Tennekes  and  Lumley  (1972)  develop  and  discuss  this  concept  in  detail. 

In  lakes,  turbulence  can  be  generated  directly  at  a  boundary 
(external  process)  or  internally  through  a  shear  instability.  The  only 
way  a  fluid  element  outside  a  turbulent  region  can  become  turbulent 
(i.e.,  acquire  vorticity)  is  through  viscous  diffusion  of  vorticity. 
Once  a  fluid  element  is  turbulent  (i.e.,  possesses  vorticity),  its 
turbulence  or  vorticity  can  be  amplified  by  the  straining  set  up  by 
neighboring  turbulence.  Because  of  this  straining,  the  vorticity  and 
energy  of  smaller  eddies  can  increase  at  the  expense  of  the  energy  of 
larger  eddies.  Energy,  therefore,  flows  from  large  to  small  eddies. 
Without  a  continual  source  of  energy  at  the  larger  scales,  turbulence 
could  not  be  maintained  and  would  rapidly  decay.  Generally,  the  larger 
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eddies  move  slowly  and  last  longer  than  the  smaller  eddies.  The  larger 
eddies  are  also  responsible  for  most  of  the  transport  of  momentum  and 
mass.  At  the  very  small  scales,  viscosity  smooths  the  velocity  fluctua¬ 
tions  by  dissipating  small-scale  motion  into  heat. 

Density  stratification  inhibits  turbulence  and  mixing.  The  reason 
for  the  decreased  mixing  in  a  zone  of  stable  density  stratification  is 
that  turbulent  eddies  lose  energy  not  only  through  viscous  dissipation, 
as  previously  described,  but  also  through  the  performance  of  work 
against  gravity  in  the  density  gradient  (Section  2.3).  The  Richardson 
number 


Ri 


(7) 


compares  the  energy  required  to  do  work  against  the  density  gradient 
(dp/dz)  with  the  energy  supply  of  the  shearing  form  (dU/dz). 

The  scattering  of  particles  by  turbulent  motion  can  be  considered 
analogous  to  molecular  diffusion  with  a  larger  turbulent  diffusion  coef 
ficient  replacing  the  molecular  diffusion  coefficient.  The  turbulence 
flux  qt  in  the  x  direction  if  therefore  given  by 


qt 


3C 

8x 


(8) 


where  Kfc 


turbulent  diffusion  coefficient, 


m  /sec. 


The  turbulent  diffusion  coefficient  is  equivalent  to  the  product  of  the 

2  1/2 

Lagrangian  length  scale  1  and  the  intensity  of  turbulence  <u  > 

Lt 


Kt  -  lt 


(9) 
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where 

<  >  *»  ensemble  mean  (l.e.f  mean  over  many  trials) 

2  1/2 

<u  >  *  root  mean  square  velocity,  m/sec 

The  Lagrangian  length  scale  is  that  distance  a  particle  must  travel 

before  it  forgets  its  initial  velocity.  This  turbulent  diffusion 

approach  is  not  valid  for  lengths  smaller  than  1^  or  times  less  than 

the  Lagrangian  time  scale  (i.e.,  the  time  required  for  a  particle  to 

move  1T  ) .  Individual  clouds  of  dimension  less  than  1T  grow  at  a 

rate  that  increases  with  size  (i.e.,  Richardson  diffusion)  and  varies 

from  cloud  to  cloud.  There  is  an  intermediate  range  where  clouds  grow 

4/3 

in  proportion  to  the  4/3  law  (i.e.,  proportional  to  1L  ),  but  this 
requires  homogeneous  turbulence  and  no  boundary  effects. 

2.4.5  Entrainment 

When  a  nonturbulent  body  of  water  is  stirred  at  a  boundary,  the 
turbulence  generated  at  that  boundary  will  advance  into  the  nonturbulent 
regions  of  the  fluid.  Entrainment  is  the  term  used  to  describe  this 
one-way  advective  type  transport  which  is  characteristic  of  most  free 
turbulent  shear  flows.  If  the  fluid  is  homogeneous,  entrainment  will 
proceed  unhindered,  and  the  thickness  of  the  layer  being  stirred  will 
increase  linearly  with  time  until  another  boundary  is  reached  or  another 
force  becomes  limiting  (e.g.,  Coriolis  effect). 

Visual  observations  in  laboratory  experiments  (Turner  1973)  showed 
the  interface  between  the  turbulent  and  nonturbulent  regions  to  be  sharp 
but  convoluted.  In  the  stratified  fluids  where  the  densities  of  the 
stirred  layer  and  nonturbulent  regions  differ,  turbulence  and  convolu¬ 
tions  at  the  interface  are  generally  suppressed  by  buoyancy  forces. 

Under  these  conditions,  entrainment  is  the  result  of  wisps  or  thin 
sheets  of  fluid  being  scoured  into  the  turbulent  region  by  turbulent 
eddies , 

The  rate  at  which  the  mean  position  of  the  interface  advances  into 
the  nonturbulent  fluid  is  usually  expressed  in  terms  of  an  entrainment 
velocity  w  .  The  entrainment  velocity  is  usually  assumed  to  be  a 


function  of  the  overall  Richardson  number 


S  Ap  h 

Ri*  -  2  (10) 

p  w* 


where 

2 

g  *  acceleration  due  to  gravity,  m/sec 

Ap  *  density  difference  between  turbulent  layer  (i.e,, 
layer  being  stirred)  and  nonturbulent  layer  (i.e., 

3 

layer  being  entrained) ,  kg/m 

h  *  depth  of  turbulent  layer,  m 
3 

p  =  density,  kg/m 

w^  *  shear  velocity  in  turbulent  layer,  m/sec 
such  that 


~  (ID 

w* 

There  is  little  reason  to  expect  Equation  11  to  be  a  simple  func¬ 
tion  because  of  the  complex  and  different  physical  processes  involved  in 
entrainment.  It  is,  however,  appealing  to  assume 


c  m;1 


(12) 


where  c  *  proportionality  constant  ,  because  the  rate  of  increase  of 
potential  energy  is  then  equal  to  the  rate  of  mechanical  energy  input. 
This  law,  however,  falls  apart  as  Ri^  -►  0  and  as  Ri^  •+■  « 

(Phillips  1977). 


2.4.6  Convection 

In  order  to  distinguish  vertical  transport  induced  by  density 
instabilities  from  advective  transport  generated  by  other  forces  or 
sources  of  energy,  convection  is  defined  as  a  buoyancy-induced  flow  that 
occurs  in  a  fluid  when  it  becomes  unstable  due  to  density  (temperature) 


differences.  If  the  mixed-layer  depth  increases  as  he,  t  is  removed  but 
no  changes  occur  in  the  thermal  profiles  below  the  mixf  d  layer,  the 
increase  in  thickness  of  the  mixed  layer  is  termed  encroachment  by 
Tennekes  and  Driedonks  (1980)  (Figure  23a).  If  instead,  some  of  the  TKE 
generated  by  cooling  at  the  air/water  interface  is  used  to  entrain  sta¬ 
bly  stratified  water  into  the  mixed  layer  and  thereby  modify  the  thermal 
profiles  below  the  mixed  layer,  the  convective  entrainment  is  termed 
penetrative  convection  (Figure  23b).  Reviews  of  penetrative  convection 
can  be  found  in  Turner  (1973)  and  Denton  (1978). 

In  order  to  explain  the  process  of  convection,  Denton  (1978) 
divided  the  water  column  into  three  zones  (Figure  24) .  The  thin  layer 
at  the  water  surface  is  the  buoyancy  production  zone.  This  layer  is  a 
molecular  diffusion  boundary  layer  that  builds  up  due  to  cooling  at  the 
air/water  interface.  At  some  critical  point,  the  layer  becomes  unstable 
and  breaks  up.  At  the  time  of  breakup,  Denton  (1978)  observed  the  for¬ 
mation  of  long,  irregular  rolls  that  eventually  broke  into  single  clumps 
of  fluid  (thermals).  Since  these  thermals  are  heavier  than  the  ambient 
fluid,  they  sink  into  the  mixed  layer.  After  the  breakup  of  the  bound¬ 
ary  layer,  a  new  layer  begins  to  form  which,  in  turn,  will  eventually 
become  unstable  and  break  up.  In  general,  the  thickness  of  the  produc¬ 
tion  layer  is  negligible  compared  to  the  total  thickness  of  the  mixed 
layer. 

The  mixed  layer,  as  implied  by  its  name,  is  uniformly  mixed  with 
respect  to  all  properties.  Dye  observations  and  continuous  temperature 
measurements  in  the  mixed  layer  indicate  that  the  thermals  generated  in 
the  production  layer  pass  through  the  mixed  layer  as  discrete  elements 
(not  plumes  as  sometimes  described),  the  net  effect  being  to  keep  the 
layer  turbulent. 

Below  the  mixed  layer  is  the  region  of  the  stably  stratified 
fluid.  The  density  gradient  is  usually  assumed  to  be  discontinuous  at 
the  interface  to  be  consistent  with  entrainment  experiments  that  report 
a  sharpening  of  the  gradient.  When  the  thermals  encounter  the  inter¬ 
face,  they  penetrate  a  short  distance  and  sometimes  rebound  (Linden 
1973).  Through  various  mechanisms  that  are  not  completely  understood, 


NATURAL  CONVECTION 


NONPENETRATIVE  PENETRATIVE 

a.  Encroachment  b.  Penetrative 

convection 

Figure  23.  Comparison  of  encroachment  and  penetrative  convection 

some  of  the  stratified  fluid  is  entrained  into  the  mixed  layer. 

According  to  Tennekeo  and  Driedonks  (1980),  10  to  50  percent  of  the 
potential  energy  released  during  cooling  is  converted  to  TKE  and  used 
for  entrainment.  Denton  (1978)  concluded  this  fraction  is  not  constant 
but  attains  a  maximum  at  a  Richardson  number  of  0.91  and  decreases  with 
increasing  and  decreasing  Richardson  number.  The  mixing  and  entrainment 
resulting  from  penetrative  convection  is  restricted  to  the  diffusive 
layer  (Figure  24),  which  has  a  thickness  of  approximately  0.2  times  the 
mixed  layer  depth. 

2.4,7  Dispersion 

Dispersion  includes  the  combined  effects  of  shear  and  diffusion. 
Shear  causes  particles  to  move  at  different  rates.  Diffusion  causes 
particles  to  move  across  the  velocity  gradient.  After  a  certain  period 
of  time,  a  particle  experiences  all  velocities  and  the  longitudinal  dis¬ 
tribution  appears  Gaussian  (Figure  25).  The  mixing  can  then  be  de¬ 
scribed  by 
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Figure  24.  Schematic  of  convective  zones 


n  ~  £C 
qd  Kd  3x 

where 

2 

K,  “  dispersion  coefficient,  m  /sec 

d  2 
=  flux  due  to  dispersion,  kg/(m  -sec) 

In  riverine  systems,  the  concentration  is  sometimes  skewed 
(Figure  25)  because  of  edge  effects.  Material  is  trapped  along  the 
banks  and  released  at  a  later  time,  causing  the  longitudinal  distri¬ 
bution  to  be  skewed  in  an  upstream  direction.  Fischer  et  al.  (1979) 
discuss  dispersion  in  detail. 
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Figure  25.  Longitudinal  dispersion 


SECTION  3.  RESERVOIR  MIXING  AND  WATER  QUALITY 


3.1  Introduction 


It  was  shown  in  Section  2.3  that  energy  is  required  to  mix  a  sta¬ 
bly  stratified  fluid.  The  energy  is  used  to  perform  work  and  raise  the 
center  of  mass.  The  major  sources  of  energy  potentially  available  to 
mix  a  reservoir  are  atmospheric  heating  (cooling),  wind,  inflows,  out¬ 
flows,  barometric  pressure,  and  gravity.  Although  several  reviews  have 
been  published  on  lake  mixing  and  hydrodynamics,  these  reviews  tend  to 
emphasize  specific  types  of  lakes  or  specific  mixing  processes.  In  gen¬ 
eral,  these  reviews  do  not  consider  all  aspects  of  mixing  and  reservoir 
operation  and  their  impacts  on  reservoir  water  quality.  For  example, 
Boyce  (1974),  Mortimer  (1974),  and  Csanady  (1975)  reviewed  the  hydro¬ 
dynamics  of  large  lakes;  Smith  (1979)  reviewed  the  hydraulics  of  iso¬ 
thermal  lakes;  Ottesen  Hansen  (1978),  Fischer  et  al.  (1979),  and 
Imberger  and  Hamblin  (1982)  reviewed  mixing  mechanisms;  Ford  and  Johnson 
(1983)  reviewed  inflow  dynamics;  Roberts  and  Dortch  (1985)  reviewed 
entrainment  of  pumped-storage  inflow  jets;  and  Imberger  (1980)  reviewed 
outflow  dynamics  and  selective  withdrawal.  In  addition,  oceanographic 
reviews  on  internal  waves  by  Garrett  and  Munk  (1979),  mixed  layer 
dynamics  by  Kraus  (1977),  and  microstructure  by  Gregg  and  Briscoe  (1979) 
are  directly  applicable  to  mixing  in  reservoirs. 

In  this  section,  reservoir  mixing  will  be  reviewed  and  discussed 
with  respect  to  energy  inputs  and  water  quality  impacts.  The  section 
concludes  with  a  summary  of  mixing  processes  that  must  be  considered 
when  developing  a  mathematical  algorithm  to  predict  water  quality 
changes . 


3.2  Description  of  Reservoir  Mixing  Processes 

3.2.1  Atmospheric  Heating  (Cooling) 

The  major  surface  heat  exchange  processes  acting  on  a  water  body 
are  summarized  in  Figure  26.  The  numbers  in  parentheses  indicate  the 
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H,  ■  SHORTWAVE  SOLAR  RADIATION  (SO  TO  400  W/m* ) 


Figure  26.  Major  heat  exchange  processes 
(after  Edinger,  Brady,  and  Geyer  1974) 


relative  magnitudes  and  ranges  of  mean  daily  values  at  midlatitude.  A 
complete  description  of  these  processes  and  their  measurement  and/or 
computation  can  be  found  in  Tennessee  Valley  Authority  (TVA)  (1972); 
Edinger,  Brady,  and  Geyer  (1974);  and  other  references.  The  relative 
magnitude  of  these  processes  depends  primarily  on  the  time  of  day  and 
year,  location  on  the  earth’s  surface  (i.e.,  latitude,  longitude,  and 
elevation),  surrounding  terrain  and  horizon  angle,  amount  of  water  in 
the  atmosphere,  and  meteorological  conditions  such  as  air  temperature, 
cloud  cover,  wind  speed  and  direction,  and  dew  point  temperature. 

Atmospheric  heating  and  cooling  impact  the  mixing  regime  in  reser¬ 
voirs  by:  (a)  adding  and  removing  heat,  thereby  changing  the  water 
density  and  RPE,  and  (b)  producing  TKE  during  periods  of  heat  loss  and 
convective  mixing.  The  seasonal  variation  in  atmospheric  heating  for 
central  Arkansas  is  shown  in  Figure  27.  It  indicates  that  periods  of 
cooling  and  convective  mixing  typically  occur  during  late  summer,  fall, 
and  winter  when  heat  losses  exceed  heat  gains.  Convective  mixing  can 
also  occur  during  the  other  periods  because  of  diurnal  heating  and  cool¬ 
ing  (Figure  28).  Hirshburg,  Goodling,  and  Maples  (1976)  analyzed  this 
mixing  process  and  developed  a  seasonal  temperature  model  based  on  this 
concept.  They  showed  that  convective  mixing  due  to  diurnal  cooling  is 
an  important  mixing  mechanism,  but  a  short  time  increment  (-1  hr)  n 
be  used  to  model  convective  mixing. 
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Figure  27.  Seasonal  variation  in  atmospheric 
heating,  central  Arkansas 

3.2.2  Wind  Processes 

The  wind  is  the  major  source  of  energy  for  many  physical  phenomena 
which  either  directly  or  indirectly  cause  mixing.  These  phenomena 
include  surface  waves,  circulation  currents,  seiches,  internal  waves, 
turbulence,  and  Langmuir  circulation. 

Wind  shear.  When  the  wind  blows  across  a  water  surface,  it 
creates  a  shear  stress  that  transmits  energy  to  the  water  body 
(Figure  29).  Part  of  this  energy  is  used  for  surface  waves,  part  for 
circulation  currents,  and  part  for  the  direct  production  of  turbulence. 
The  magnitude  of  the  shear  stress  is  determined  from 
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(14) 


where 

2 

i  =  surface  shear  stress,  kg/(m-sec  ) 

s  3 

p  =  density  of  air,  1.177  kg/m 

3  -3 

*=  drag  coefficient  (order  of  10  )  ,  dimensionless 

W  =  wind  speed  at  a  specified  elevation,  m/sec. 
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Figure  28,  Diurnal  heating  and  cooling  cycles 
for  central  Arkansas  (Continued) 
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Figure  28.  (Concluded) 
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Figure  29.  Wind  shear  on  water  surface 
(after  Ottesen  Hansen  1978) 

The  drag  coefficient  can  be  determined  from  empirical  formulas  developed 
by  Safaie  (1978),  Wu  (1980),  and  others.  The  rate  of  energy  input  to 
the  water  body  is 
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d(TKE) 

dt 


(15) 


t  dA 
s 


where 

A  =  water  surface  area,  m^ 
s 

ug  =  surface  drift  velocity,  m/sec 

The  surface  drift  velocity  is  usually  assumed  to  scale  with  the  wind 
friction  velocity  in  water  (w^) ,  which  is  defined  by 
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u 
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(16) 


3 

where  p  =  water  density,  kg/m  . 
w 

From  Equations  14-16,  it  is  evident  that  the  energy  available  to 
create  waves,  turbulence,  and  currents  and,  hence,  mixing  is  dependent 
on  the  wind  speed  to  the  third  power  and  on  the  surface  $rea  of  the 
reservoir.  Depending  on  the  surrounding  terrain  and  shape  of  the  reser¬ 
voir,  measured  wind  speeds  and  surface  areas  may  not  be  appropriate  for 
use  in  Equations  14  and  15.  If  the  terrain  surrounding  a  reservoir 
consists  of  high  hills,  trees,  bluffs,  etc.,  the  terrain  may  shelter  the 
water  surface  from  the  wind  force.  As  a  rule  of  thumb,  the  sheltering 
effects  extend  into  the  lake  a  distance  of  approximately  eight  times  the 
vertical  elevation  (Ford  1976).  The  wind  speed  can  also  increase  over  a 
water  surface  because  of  the  decrease  in  surface  roughness  going  from 
land  to  water.  Using  boundary  layer  theory.  Ford  and  Stefan  (1980b) 
developed  a  relationship  to  predict  the  increase  in  wind  speed  based  on 
fetch  and  roughness  lengths.  For  complete  reviews  on  this  subject,  the 
reader  is  referred  to  Haugen  (1973)  and  Smith  (1975).  Dendritic  reser¬ 
voirs  with  many  islands  may  experience  funneling  and  variability  in  wind 
speed  and  direction  which  further  complicate  the  determination  of  wind 
shear . 
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Surface  waves.  The  significance  of  wave  action  with  respect  to 
mixing  results  from  the  orbital  motion  of  water  particles  beneath  the 
water  surface  (Figure  30) .  In  isothermal  lakes  or  in  the  mixed  layer  of 
a  stratified  lake,  the  depth  of  water  affected  by  wave  action  is  approx¬ 
imately  one-half  the  wave  length  (X)  (Smith  1979).  This  depth  is  also 
used  to  define  the  shore  zone  or  that  region  of  a  lake  where  the  wave 
form  is  distorted  by  bottom  interference. 
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Figure  30.  Schematic  of  surface  waves  (after  Smith  1979) 


The  waves  become  unstable  and  break  when  the  water  depth  is  equal  to 
four-thirds  the  wave  height  (H) .  This  point  is  the  break  line  and 
defines  the  swash  zone  as  the  area  where  the  waves  break  out  of  shore. 
Since  the  orbitals  are  not  truly  circular  and  do  not  completely  close, 
there  is  a  net  transport  by  wave  motion  in  the  direction  of  the  wind. 
This  transport  is  usually  negligible  compared  to  the  transport  by  cur¬ 
rents  (Smith  1979). 

3.2.3  Circulation  Currents 

General.  Circulation  currents  are  averaged  water  movements  con¬ 
trolled  by  Coriolis,  external,  and  friction  forces.  Circulation  pat¬ 
terns  In  reservoirs  are  complicated  by  basin  configurations,  density 
stratification,  fluctuating  wind  speeds  and  directions,  reservoir 
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operating  plans  (e.g.,  hydropower  and  flood  control),  and  secondary 
circulations  from  differential  heating  and  convective  currents.  In 
addition,  the  principle  of  continuity  requires  that  any  movement  of 
water  out  of  any  part  of  a  reservoir  be  balanced  by  the  return  of  an 
equal  volume  of  water  or  a  change  in  water  surface  level. 

Coriolis  forces,  which  result  from  the  earth's  rotation  about  its 
axis,  cause  currents  in  the  northern  hemisphere  to  deflect  to  the  right 
when  looking  in  the  direction  of  the  flow.  For  example,  Pharo  and 
Carmack  (1979)  determined  that  Coriolis  effects  caused  the  Thompson 
River  inflow  to  Kamloops  Lake,  British  Columbia,  to  be  deflected  along 
the  shore.  According  to  Mortimer  (1974),  Coriolis  forces  or  rotational 
effects  become  significant  when  the  width  of  a  lake  is  greater  than  5r 
and  dominant  when  the  lake  width  is  greater  than  20r.  The  radius  of  the 
inertial  circle  r  is  defined  by 


r 


(17) 


where 

U  *  water  velocity,  m/sec 
f  *  Coriolis  parameter  (f  -  2w  sin  0) 

}  =  angular  velocity  of  the  earth’s  rotation 
(7.29  x  lcf5  rad/sec) 

0  *  latitude 

For  a  typical  velocity  U  of  0.1  m/sec  and  a  latitude  0  of 
40  deg,  r  1.06  km  and  rotational  effects  become  significant  for  lakes 
with  widths  greater  than  5  km.  Because  r  is  inversely  related  to  the 
sine  of  the  latitude  0  ,  rotational  effects  become  significant  at  a 
smaller  width  with  increasing  latitude.  In  large  lakes,  which  are 
dominated  by  rotational  effects,  the  classical  Ekman  spiral  can  limit 
the  depth  of  the  upper  mixed  layer  to  the  Ekman  layer  depth  and  thereby 
affect  the  mixing  regime  (Smith  1979)  . 

The  wind  is  a  major  external  force  governing  circulation  patterns 
in  lakes  and  reservoirs.  As  previously  discussed,  the  fraction  of  the 
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wind  shear  that  is  used  to  generate  currents  is  unknown.  It  is,  how¬ 
ever,  usually  assumed  that  the  surface  current  speed  is  1  to  3  percent 
of  the  wind  speed.  Bengtsson  (1978)  and  others  have  found  that: 

(a)  the  percentage  is  not  a  constant  even  for  a  specific  case,  (b)  the 
percentage  decreases  with  increasing  wind  speed,  and  (c)  the  percentage 
increases  with  increasing  lake  dimensions.  Figure  31  illustrates  the 
complex  surface  current  patterns  that  can  develop  in  a  lake.  Although 


SCALE  km 

Figure  31.  Surface  currents  in  Lake  Ivo, 
Sweden  (after  Bengtsson  1978) 


it  Is  usually  assumed  that  the  current  speeds  decline  exponentially  with 
depth,  the  actual  shape  of  the  current  vertical  profile  will  depend  on 
the  nature  of  all  external  forces,  stratification,  inflow  and  outflow 
distributions,  and  the  location  of  all  lake  boundaries.  Friction  at  the 
lake  boundaries  causes  a  boundary  layer  to  form  that  can  be  described  by 


a  logarithmic  velocity  profile  such  that  the  current  speed  at  the  bed  is 
zero.  The  random  variability  of  the  wind  and  the  geometrical  complexi¬ 
ties  of  reservoir  basins  result  in  a  temporally  changing  and  spatially 
nonuniform  water.  Smith  (1979)  discusses  many  of  these  and  other 
factors  controlling  lake  circulation  in  detail. 

In  reservoirs,  the  existence  of  well-defined  inflow  density  cur¬ 
rents  (Ford  and  Johnson  1983),  withdrawal  zones  (Bohan  and  Grace  1973, 
Xmberger  1980),  and  possibly  pumpback  jets  (Roberts  1981)  probably  sig¬ 
nificantly  perturbs  classical  wind-generated  circulation  patterns,  but 
little  is  known  of  these  interactions.  In  addition,  reservoir  operating 
procedures  that  result  in  pool-level  fluctuations  also  generate  currents 
that  must  interact  with  the  wind-generated  current  field.  Again,  little 
is  known  about  the  synergistic  effects  of  these  currents  and  mixing 
mechanisms  except  that,  once  water  is  set  in  motion,  it  cannot  respond 
instantaneously  to  a  change  in  flow  regime. 

Seiches .  When  a  steady  wind  blows  across  a  water  surface,  water 
is  transported  to  the  windward  side  of  the  lake  and  the  water  surface 
slope  is  determined  by 
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where 

S  =  water  surface  slope 
AD  =  setup  of  water  surface,  L 
L  =  length  of  lake,  L 

2 

t  *  surface  shear  stress,  F/L 
s 

D  =  water  depth,  L 


In  this  equation,  the  wind  force  is  balanced  by  the  change  in  hydro¬ 
static  pressure.  If  the  water  body  is  stratified  such  that  it  can  be 
represented  as  a  two-layered  system,  the  interface  responds  to  the  wind 
force  by  tilting  in  the  opposite  direction  (Figure  32)  with  a  slope  of 
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(19) 
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where 


S.  *  interface  slope 
■1 

Ap  ■  density  difference  between  two  layers,  m/L' 
*  thickness  of  top  layer,  L 


AD/2 


If  the  wind  force  (t  )  is  sufficiently  large,  the  interface  slope  may 
become  large  enough  so  that  the  lower  layer  is  exposed  to  the  water 
surface.  This  phenomenon  is  termed  upwelling  and  can  result  in  the 
hypolimnetic  water  being  mixed  into  the  surface  waters  (Blanton  1973). 

When  the  wind  force  is  removed,  the  potential  energy  associated 
with  the  water  surface  and  Interface  is  converted  into  kinetic  energy 
and  flow.  The  result  is  an  oscillating  motion  called  seiching,  which 
decays  with  time.  The  time  for  one  cycle  of  the  surface  seiche  (t.e., 
the  seiche  period)  is 
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where  the  terms  are  as  previously  defined.  A  seiche  Is  one  Indirect 
mechanism  for  energy  from  the  wind  to  become  available  for  mixing  in  the 
hypolimnion  of  a  lake.  A  more  detailed  review  of  seiches  and  upwelli  , 
is  found  in  Monismlth  (1983). 

3.2.4  Internal  Waves 

All  stably  stratified  fluids  possess  the  ability  to  sustain 
internal  wave  motion.  The  waves  can  be  generated  locally  by  turbulence 
or  externally  by  an  outside  perturbation  such  as  hydropower  operation. 
Internal  waves  transport  momentum  and  can  exist  without  breaking  and 
forming  turbulence.  It  is  only  after  they  break  that  turbulence  is 
generated  and  mixing  is  possible. 

Internal  waves  Impact  the  mixing  regime  because  they  radiate 
energy  away  from  the  place  where  they  were  generated.  They  can  reduce 
entrainment  into  the  mixed  layer  by  creating  a  local  energy  loss  from 
the  mixed  layer  (Kantha,  Phillips,  and  Azad  1977).  Internal  waves  can 
increase  mixing  by  transporting  energy  into  a  region  where  it  would  not 
otherwise  be  available. 

Turbulence.  Part  of  the  energy  input  from  the  wind  shear  is  used 
to  directly  produce  turbulence  and  TKE.  This  turbulence  interacts  with 
turbulence  produced  by  breaking  waves  (surface  and  internal) ,  cooling  at 
the  air/water  interface,  and  current  shear  to  keep  the  upper  layer 
(epilimnion)  well  mixed  and  to  entrain  metalimnetic  water  into  the 
epilimnicn.  If  the  turbulence  is  generated  at  a  length  scale  less  than 
the  thickness  of  the  mixed  layer,  viscous  dissipation  with  depth  will 
result  in  a  negligible  fraction  of  the  energy  being  available  for 
entrainment  at  the  epilimnetic/metalimnetic  interface.  If  there  is  a 
simultaneous  heat  (buoyancy)  flux  across  the  air/water  interface,  all  of 
the  wind-generated  TKE  may  be  used  to  maintain  the  well-mixed  layer  and 
not  be  available  for  additional  entrainment. 

Langmuir  circulation.  Langmuir  circulations  are  counter-rotating 
concentric  vortices  in  the  surface  layers  of  water  bodies  (Figure  33). 
Their  existence  is  readily  observed  because  materials  accumulating  in 
the  zones  of  convergence  appear  as  streaks  or  bands  on  the  water 
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Figure  33.  Schematic  of  a  Langmuir  circulation 
cell  (after  Pollard  1977) 


surface.  In  some  instances,  compressed  films  dampen  the  capillary 
waves,  giving  the  streaks  a  smoother  appearance.  Langmuir  (1938)  was 
the  first  to  show  that  the  commonly  observed  surface  streaks  or  windrows 
were  a  manifestation  of  a  more  complex  circulation  pattern  related  to 
the  wind.  Langmuir  also  concluded  that  the  vortices  were  responsible 
for  the  formation  of  stratification  and  the  maintenance  of  a  well-mixed 
layer.  Recent  laboratory  and  theoretical  studies  have  shown  that  both 
waves  and  current  shear  are  required  for  the  Langmuir  cells  to  form 
(Leibovich  1983).  Earlier  models  that  related  the  formation  of  Langmuir 
circulation  to  instabilities  and  thermal  convection  are  now  considered 
invalid.  In  a  comprehensive  review  of  Langmuir  circulation,  Leibovich 
(1983)  summarizes  the  results  from  several  field  studies.  These 
*  include: 

i 
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a.  Langmuir  cells  form  and  align  within  a  few  degrees  of  the  wind 
direction, 

b.  Cells  form  within  a  few  minutes  of  onset  of  a  wind  of  3  m/sec 
or  faster, 

c.  Spacing  of  windrows  varies  from  1  m  to  hundreds  of  metres. 

d.  Depth  of  penetration  of  cells  is  limited  to  the  first 
significant  density  gradient. 

e.  The  spacing  between  the  largest  windrows  is  approximately 
twice  the  depth  of  the  cells. 


Although  there  is  evidence  to  suggest  Langmuir  cells  are  important 
mixing  mechanisms  in  reservoirs,  there  is  no  direct  proof  (Leibovich 
1983). 


3.2.5  Inflows 

Since  tributary  inflow  density  usually  differs  from  the  density  of 
the  reservoir  surface  waters,  inflows  enter  and  move  through  reservoirs 
as  density  currents.  Depending  on  the  sign  and  magnitude  of  this  den¬ 
sity  difference,  density  currents  can  enter  the  epilimnion,  metalimnion, 
or  hypolimnion  (Figure  34)  and  thereby  modify  the  mixing  regime  in  these 
regions  since  inflows  are  a  source  of  both  buoyant  potential  energy  and 
turbulent  kinetic  energy.  The  rate  of  the  PE  input  is  dependent  on  the 
inflow  density,  the  flow  rate,  and  the  in-lake  density  distribution: 

-  It  (m8H)  ■  Vzj  <21) 


where 

Ap .  ■  density  difference  between  inflow  and  reservoir  surface 
3  3 

waters,  kg/m 

3 

Qj  =  inflow  rate,  m  /sec 

Z.  =  difference  in  depth  between  center  of  mass  of  inflow  and 
^  center  of  mass  of  inflow  placement,  m 
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Figure  34.  Inflow  density  currents 
The  rate  of  TKE  input  is  also  dependent  on  the  flow  rate: 

^  (IKE)  -  k  ({  mU2)  -  i.jQj  U2 


where 


inflow  velocity  *  Q./A  ,  m/sec 
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inflow  density,  kg/ni 
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cross-sectional  area  perpendicular  to  inflow,  m 
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Since  the  mixing  resulting  from  inflows  is  dependent  on  the  energy 
input,  the  magnitude  of  inflow  mixing  is  highly  dependent  on  the  flow 
rate.  A  complete  analysis  of  reservoir  inflows  is  found  in  Ford  and 
Johnson  (1983). 

When  the  inflow  density  is  less  than  the  surface  water  density, 
the  inflow  will  float  on  the  water  surface  as  an  overflow.  This  con¬ 
dition  typically  occurs  in  the  spring  when  inflows  warm  more  rapidly 
than  reservoirs.  In  an  overflow,  the  excess  hydrostatic  pressure  in  the 
density  current  causes  the  current  to  flow  in  all  directions  not 
obstructed  by  boundaries.  Overflows  are  susceptible  to  mixing  from 
wind-induced  mechanisms  and  diurnal  heating  and  cooling  processes.  Wind 
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shear  can  direct  the  overflow  into  a  cove  or  prevent  it  from  moving 
downstream.  Vertical  mixing  in  the  mixed  layer  can  be  enhanced  with  an 
overflow  if  the  density  difference  is  small  and  flow  rate  is  large 
(i.e.,  large  TKE  input)  or  suppressed  if  the  flow  rate  is  small  and 
density  difference  is  large  (this  situation  frequently  exists  with 
thermal  discharges) . 

If  the  inflow  density  is  greater  than  the  water  surface  density, 
the  inflow  will  push  the  reservoir  water  ahead  until  the  buoyancy  forces 
dominate  and  the  inflow  plunges  beneath  the  water  surface.  The  plunge 
point  Is  sometimes  made  visible  because  of  turbidity  or  the  accumulation 
of  floating  debris,  which  may  indicate  a  stagnation  point.  The  location 
of  the  plunge  point  is  determined  by  a  balance  between  the  stream  momen¬ 
tum,  the  pressure  gradient  across  the  interface  separating  the  river  and 
reservoir  waters,  and  the  resisting  shear  forces.  Some  mixing  (termed 
initial  mixing)  occurs  at  the  plunge  point  because  of  the  large  eddies 
formed  by  flow  reversals  and  pooling  of  the  inflowing  water  (Akiyama  and 
Stefan  1981).  Ford  and  Johnson  (1983)  estimated  this  mixing  to  be  on 
the  order  of  25  percent  of  the  inflow  volume.  Knapp  (1942)  noticed  that 
the  flow  in  the  vicinity  of  the  plunge  point  occurred  at  the  bottom  of 
this  pooled  mixing  zone  (Figure  35).  Ford,  Johnson,  and  Monismith 
(1980)  and  Kennedy,  Gunkel,  and  Carlile  (1983)  substantiated  this  pool¬ 
ing  phenomenon  during  dye  studies  on  DeGray  Lake,  Arkansas,  and  West 
Point  Lake,  Georgia,  respectively,  when  the  dye  clouds  appeared  to  have 
stalled  at  the  plunge  point.  The  location  of  the  plunge  point  can  also 
be  influenced  by  morphological  factors.  Changes  in  the  bed  slope  (e.g., 
due  to  sediment  deposition),  bed  friction,  and  cross-sectional  area  may 


Figure  35.  Flow  in  the  vicinity  of 
the  plunge  point 
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affect  location  of  the  plunge  point.  For  a  river  entering  a  wide  lake, 
the  plunge  point  may  actually  be  a  point  rather  than  a  line. 

After  the  inflow  plunges,  it  follows  the  old  river  channel 
(thalweg)  as  an  underflow.  The  speed  and  thickness  of  the  underflow  is 
determined  by  a  flow  balance  between  the  shear  forces  and  the  accelera¬ 
tion  due  to  gravity  (i.e.,  gradually  varying  flow  theory).  An  underflow 
will  entrain  overlying  reservoir  water  due  to  shear  and  turbulence  gen¬ 
erated  by  bottom  roughness.  Changes  in  the  underflow  density  from 
entrainment  must  be  quantified  before  a  density  interflow  or  intrusion 
can  be  analyzed,  or  estimates  of  the  vertical  placement  of  an  interflow 
can  be  incorrect. 

An  interflow  or  intrusion  occurs  when  a  density  current  leaves  the 
river  bottom  and  propagates  horizontally  into  a  stratified  body  of 
water.  Intrusions  differ  from  overflows  and  underflows  because  an 
intrusion  moves  through  a  reservoir  at  a  elevation  where  the  intrusion 
and  reservoir  densities  are  similar.  Intrusions  require  a  continuous 
inflow  and/or  outflow  for  movement,  or  they  stall  and  collapse  (i.e., 
dissipate).  Turbulence  is  usually  quickly  dissipated  in  an  intrusion 
since  the  metalimnetic  density  gradient  creates  strong  buoyancy  forces 
which  inhibit  mixing.  Mixing  still  occurs,  however,  because  of  the  flow 
gradient . 

3.2.6  Outflows 

Manmade  reservoirs  differ  from  natural  lakes  in  many  respects,  but 
especially  with  respect  to  the  importance  of  outflows.  In  natural 
lakes,  outflows  are  usually  uncontrolled  surface  withdrawals  with  the 
outflow  rate  dependent  on  the  water  surface  elevation.  In  reservoirs, 
outflows  are  usually  regulated  by  gates  and/or  a  structure,  and  the 
withdrawal  depth  is  not  necessarily  near  the  water  surface.  For  exam¬ 
ple,  many  flood  control  projects  have  bottom  outlets  while  multipurpose 
projects  may  have  multilevel  outlets  (i.e.,  selective  withdrawal 
structures . 

Outflows  contribute  to  the  in-lake  mixing  regime  through  with¬ 
drawal  currents  and  turbulence  associated  with  them.  Withdrawal 


currents  are  generated  because  the  PE  associated  with  the  water  level 
(i.e.f  the  head)  is  converted  into  TKE  and  causes  water  to  flow  through 
an  outlet  (see  Figure  36).  The  characteristics  of  the  withdrawal  cur¬ 
rents  will  be  dependent  on  the  withdrawal  rate  (i.e.,  gate  or  structure 
setting),  ambient  in-lake  stratification,  outlet  location,  and  lake 
geometry.  If  the  outlet  is  located  in  the  hypolimnion  of  a  strongly 
stratified  lake,  the  withdrawal  zone  may  also  be  limited  to  the 
hypolimnion  of  the  lake.  Withdrawal  zones  can  also  be  limited  to  the 
epilimnion  or  metalimnion  of  a  lake  depending  on  the  outlet  location, 
flow  rate,  and  ambient  stratification.  Withdrawal  zones  can  be  computed 
using  formulas  provided  by  Bohan  and  Grace  (1973),  Imberger  (1980),  and 
Davis,  Holland,  and  Wilhelms  (1985). 


OUTFLOW 


£  (KE)  =  V2V2dJ?  s  VzV7PQ 
dt  dt 

Figure  36.  Energy  conversion  associated  with  an  outlet 

3.2.7  Barometric  Pressure  and  Gravity 

On  large  lakes,  horizontal  variations  in  barometric  pressure  can 
result  in  surface  seiching.  As  indicated  in  Section  3.2.2.,  seiching 
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can  impact  the  mixing  regime  as  internal  waves  and  seiches  form  and 
break.  Tidal  movement  resulting  from  the  gravitational  attraction  of 
the  moon  and  the  sun  has  been  observed  in  very  large  lakes  (e.g.,  Lake 
Superior,  Minnesota  and  Wisconsin,  and  Lake  Baikal,  Siberia)  (Wetzel 
1983)  but  is  probably  negligible  in  most  reservoirs. 

3.2.8  Summary 

In  the  preceding  discussion  of  reservoir  mixing,  it  was  shown  that 
the  observed  thermal  structure  in  a  reservoir  results  from  the  cumula¬ 
tive  effects  of  a  number  of  complex,  nonlinear,  interdependent  mixing 
processes.  Although  the  sources  of  TKE  for  mixing  are  limited  (i.e., 
solar  and  atmospheric  heating,  wind,  inflows,  and  outflows  (including 
outlet  location))  and  well  known,  quantitative  knowledge  of  specific 
mixing  mechanisms  is  limited.  TKE  budgets  for  the  entire  lake  do,  how¬ 
ever,  clearly  expose  the  physical  causes  of  observed  motion  and  thermal 
stratification. 

The  observed  temperature  structure  in  a  lake  is  actually  a  signa¬ 
ture  of  past  mixing  events.  The  relative  timing  of  these  events  is 
important.  For  example,  the  occurrence  of  spring  floods  after  stratifi¬ 
cation  forms  can  result  in  a  significantly  different  thermal  structure 
than  if  the  floods  occur  prior  to  stratification.  Identification  of  the 
mixing  events  is  crucial  to  understanding  the  thermal  structure  of  a 
lake . 


3.3  Influence  of  Mixing  on  Reservoir  Water  Quality 

While  the  mixing  mechanisms  of  inflows  and  outflows,  upwelling  and 
seiches,  turbulence,  and  others  are  highly  interactive,  the  influence  of 
each  of  these  processes  on  chemical  and  biological  processes  will  be 
discussed  separately  for  convenience.  It  must  be  recognized,  however, 
that  reservoir  water  quality  is  a  function  of  the  interactions  among  the 
dynamic  mixing  processes  and  of  the  chemical  and  biological  responses  to 
mixing . 
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3.3.1  Inflows  and  Outflows 


Reservoir  inflows  usually  represent  the  major  contributions  of 
both  dissolved  and  particulate  constituents  to  the  reservoir  mass 
budget.  Ground-water  and  atmospheric  contributions  are  generally 
minimal . 

As  an  inflow  enters  the  lake,  velocity  and  turbulence  decrease, 
coarse  particulate  or  suspended  materials  settle,  and  the  water  clarity 
increases.  Turbulence  generated  through  bottom  shear  generally  is  suf¬ 
ficient  to  maintain  silt-  and  clay-sized  particles  in  suspension. 
Nutrients,  bacteria,  and  organic  constituents  may  be  transported  into 
the  reservoir  since  these  constituents  generally  sorb  to  particles  in 
the  clay-silt  size  range  (Frink  1969;  Pita  and  Hyne  1974;  Sharpley  and 
Segers  1979;  De  Pinto,  Young,  and  Martin  1982).  Depending  on  the 
density  structure  of  the  lake  and  inflow,  the  inflow  can  enter  the 
epilimnion,  metalimnion,  or  hypolimnion  (Figure  34). 

An  overflow  situation  may  have  the  greatest  initial  influence  on 
reservoir  water  quality  by  introducing  oxygen-demanding  material,  nutri¬ 
ents,  and  bacteria  directly  into  the  surface  waters  or  euphotic  zone. 
Available  nutrients  can  be  assimilated  by  phytoplankton  and  may  stimu¬ 
late  plankton  blooms;  bacteria  concentrations  may  exceed  body  contact 
recreation  standards;  and  oxygen-demanding  material  may  depress  mixed- 
layer  dissolved  oxygen  (DO)  concentrations.  Even  though  the  euphotic 
zone  thickness  would  probably  be  decreased  due  to  increased  turbidity, 
circulation  and  mixing  in  the  mixed  layer  would  circulate  phytoplankton 
between  the  nutrient-rich  inflow  and  euphotic  zone.  With  overflow  con¬ 
ditions,  phytoplankton  concentration  would  probably  attain  a  maximum  in 
the  headwater  regions  and  decrease  in  the  downstream  direction  (Kimmel, 
Lind,  and  Paulson  1984). 

Interflows  and  underflows  have  similar  influences  on  water  qual¬ 
ity.  At  the  plunge  point,  the  inflow  waters  sink  to  form  a  density  cur¬ 
rent  flowing  along  the  bottom.  Bottom-generated  turbulence  results  in 
entrainment  of  surface  waters  into  the  density  current.  This  also  may 
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introduce  nutrients,  organic  matter,  and  bacteria  into  the  overlying 
surface  water  and  euphotic  zone. 

Once  the  density  current  enters  the  metalimnion  as  an  interflow, 
the  inflowing  constituents  are  temporarily  isolated  from  the  mixed 
layer.  These  constituents  may  be  entrained  into  the  mixed  layer  during 
storm  events  if  wind  gusts  deepen  the  mixed  layer  or  as  the  thermocline 
erodes  due  to  penetrative  convective  mixing.  Since  the  metalimnion  is 
isolated  from  oxygen  exchange  with  the  surface,  oxygen-demanding  mate¬ 
rial  in  the  interflow  may  depress  metalimnetic  DO  concentrations.  In 
most  stratified  reservoirs,  a  metalimnetic  DO  minimum  typically  develops 
with  anoxic  conditions  developing  in  some  systems  (Hannan  and  Cole 
1984) . 

Underflows  not  only  introduce  inflowing  constituents  into  the 
hypolimnion  but  also  increase  turbulent  diffusion  across  the  sediment/ 
water  interface  because  of  bottom  shear.  This  can  increase  oxygen 
depletion  by  increasing  oxygen  transfer  across  the  diffusion-limited 
sediment/water  interface.  It  also  can  increase  the  transfer  of  reduced 
and  resolubilized  constituents  such  as  ammonia,  phosphorus,  iron,  man¬ 
ganese,  and  hydrogen  sulfide  from  the  sediment  into  the  overlying  water 
column. 

Interflows  and  underflows  have  been  proposed  as  major  transport 
processes  for  reduced  and  resolubilized  species  such  as  manganese  and 
phosphorus  from  the  upper  reservoir  to  downstream  areas  (Nix  1981,  1984; 
Davison,  Woof,  and  Rigg  1982).  Since  anoxic  conditions  generally  begin 
in  the  upstream  area  of  the  reservoir,  inflow  mixing  processes  can 
promote  the  movement  of  this  anoxic  front  further  into  the  reservoir. 
Many  dissolved  nutrient  species  are  in  a  readily  available  form  for 
phytoplankton  assimilation. 

Underflows  also  can  improve  reservoir  water  quality.  Cold,  well- 
oxygenated  outflows  from  Lake  Ouachita,  Arkansas,  maintain  an  aerobic 
hypolimnion  in  the  old  thalweg  of  Lake  Hamilton,  Arkansas,  downstream, 
even  though  the  Lake  Hamilton  coves  become  anoxic. 

The  above  situation  can  be  altered  depending  on  how  the  project  is 
operated  (i.e.,  how  water  is  released).  If  a  reservoir  has  a  multilevel 
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withdrawal  structure,  water  may  be  released  at  a  similar  elevation  as 
the  inflow  enters  the  reservoir.  Under  such  conditions  the  inflow 
waters  can  be  short-circuited  through  the  reservoir.  If  the  residence 
time  of  a  density  overflow  is  less  than  the  response  time  (i.e.,  dou¬ 
bling  time)  of  phytoplankton  under  natural  conditions,  phytoplankton 
blooms  may  be  minimized.  Lund,  Mackereth,  and  Mortimer  (1963)  report 
natural  doubling  times  for  Asterionella  formosa  to  be  5  to  7  days. 
Doubling  times  for  other  species  under  natural  conditions  are  on  the 
order  of  5  to  30  days  (Ford  and  Thornton  1979).  Therefore,  when  the 
residence  time  or  short-circuited  residence  time  is  less  than  approxi¬ 
mately  5  days,  phytoplankton  washout  is  probable.  For  Lake  Windermere, 
England,  a  natural  lake,  Lund  (1950)  concluded  losses  due  to  overflow 
were  approximately  compensated  for  by  inflowing  nutrients.  Lund, 
Mackereth,  and  Mortimer  (1963)  found  outflow  losses  to  be  the  major  loss 
for  Asterionella.  Residence  times  of  5  days  or  less  for  overflows  imply 
that  the  inflow  water  quality  will  dominate  the  mixed  layer  quality. 

If  the  outflow  is  released  from  a  different  level  than  the  inflow 
enters,  mixing  within  the  reservoir  may  be  increased.  Transport  of 
nutrients  from  the  metalimnion  to  the  epilimnion  and  euphotic  zone,  for 
example,  can  be  increased  when  inflows  enter  as  interflows  or  underflows 
and  the  project  is  operated  with  surface  discharge.  Bottom  withdrawal 
may  promote  movement  of  interflows  and  underflows  through  the  reservoir 
and  minimize  the  interaction  between  the  inflowing  constituents  and  the 
euphotic  zone.  Bottom  withdrawal  also  may  purge  an  anoxic  hypolimnion 
of  reduced  and  resolubilized  constituents.  Dunst  (1974)  found  that  bot¬ 
tom  discharge  from  several  Wisconsin  reservoirs  with  anoxic  hypolimnia 
significantly  reduced  the  hypolimnetic  phosphorus  load.  During  fall 
overturn,  then,  this  phosphorus  load  would  not  be  available  for  mixing 
throughout  the  water  column  and  phytoplankton  assimilation. 

Many  reservoirs  are  operated  to  store  water  for  downstream  flood 
control.  The  pool  elevation  therefore  rises  to  accommodate  the  inflow¬ 
ing  water.  Water  also  flows  into  littoral  zones  and  coves,  carrying 
with  it  inflowing  constituents.  These  constituents  can  remain  within 
these  zones  and  coves  and  impact  water  quality  later  in  the  season. 
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Nutrients,  for  example,  may  be  assimilated  by  phytoplankton,  and  organic 
loadings  may  contribute  to  the  onset  of  anoxic  conditions.  When  there 
is  a  net  outflow  from  a  reservoir,  there  is  also  a  net  transport  from 
the  littoral  zones  and  coves  into  the  main  body  of  the  reservoir. 
Increased  horizontal  transport  also  results  in  increased  vertical 
transport. 

3,3.2  Upwelling,  .‘’eiches,  and  Internal  Waves 

The  water  quality  significance  of  upwelling  and  seiches  is 
severalfold.  First,  upwelling  is  one  mechanism  whereby  the  waters  of 
the  metalimnion,  and  possibly  the  hypolimnion,  are  temporarily  exposed 
to  the  surface.  Upwelling  of  anoxic  hypolimnetic  water  near  the  dam  in 
C.  J.  Brown  Lake,  Ohio,  depressed  surface  DO  concentrations  to  about 
2  mg/Jt  with  increased  surface  nutrient  concentrations  and  hydrogen 
sulfide  odors.  This  upwelling  occurred  because  of  the  passage  of  a 
large  storm  front.  The  colder,  nutrient-rich  water  may  initiate  a 
phytoplankton  bloom  of  different  composition  than  previously  occurred 
within  the  lake. 

Second,  seiches  may  alter  the  mixed-layer  depth  and  light  regime 
to  which  plankton  are  exposed.  Horizontal  variations  are  also  expected 
because  the  wind  may  concentrate  phytoplankton  at  the  downwind  end  of 
the  lake,  creating  differences  in  mixed-layer  depths  and  light  penetra¬ 
tion  within  the  mixed  layer.  The  correlation  between  upwelling  and 
phytoplankton  blooms  is  well  established  for  coastal  zones  (Fogg  1975, 
Rao  1977)  and  Is  also  found  in  natural  lakes  (Coulter  1968).  The  high 
rate  of  production  throughout  the  year  in  Lake  Victoria,  Africa,  is 
attributed  to  efficient  mixing  resulting  from  seiche  action  (Fogg  1975). 

Seiche  motion  may  indirectly  initiate  hypolimnetic  mixing.  As  the 
water  oscillates  In  the  reservoir,  bottom  shear  may  create  turbulence 
that  increases  mixing  in  the  hypolimnion.  This  mixing  may  increase 
transfer  of  DO  and  reduced  constituents  across  the  sediment/water 
interface . 

Seiche  motion  may  result  direccly  from  project  operation,  par¬ 
ticularly  in  peaking  hydropower  reservoirs.  During  generation,  water 
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movement  is  established  toward  the  dam.  When  generation  ceases,  water 
continues  to  move  toward  the  dam  and  piles  up  at  the  dam  with  the 
establishment  of  oscillations  of  this  water  as  the  water  surface  seeks 
equilibrium.  Seiche  activity  may  explain  the  observed  hypolimnetic  DO 
fluctuations  in  Lake  Sidney  Lanier,  Georgia.  Hypolimnetic  DO  concentra¬ 
tions  at  a  given  depth  near  Buford  Dam  increased  about  2  mg/i  around 
1000  hr  from  Tuesday  through  Saturday  but  were  not  evident  on  Sunday  or 
Monday.  Over  long  weekends  with  no  generation,  the  increase  was  delayed 
until  the  day  following  the  initiation  of  generation.  Apparently,  a 
seiche  was  created  through  generation  with  a  frequency  in  phase  with  the 
generation  cycle.  This  seiche  motion  created  by  hydropower  generation 
may  also  increase,  hypolimnetic  mixing  through  bottom  shear. 

While  internal  waves  do  not  promote  mixing  unless  these  waves 
break,  the  vertical  displacement  of  water  by  internal  waves  can  signifi¬ 
cantly  alter  the  light  regime  of  phytoplankton  in  the  metalimnion. 
Metalimnetic  nutrient  concentrations  may  be  higher  than  in  the  mixed 
layer  due  to  interflow  loadings  and_microbial  decomposition  of  organic 
matter  settling  from  the  epilimnion.  Light  may  be  the  factor  limiting 
phytoplankton  production  in  the  metalimnion.  Internal  waves  with  an 
amplitude  of  2  to  15  m  and  a  period  of  15  min  to  12  hr  have  been 
observed  in  coastal  areas  (Reid  1956,  LaFond  1962).  Armstrong  and 
LaFond  (1966)  have  observed  vertical  displacements  of  nutrient  layers  by 
these  waves.  As  the  metalimnion  is  displaced  by  the  internal  waves, 
this  layer  may  enter  the  ei  photic  zone  where  available  light  may  stimu¬ 
late  plankcon  production  (Denman  and  Gargett  1983).  Plankton  production 
may  increase  metalimnetic  DO  concentrations  during  the  day,  but 
increased  respiration  at  night  may  deplete  DO  concentrations  below 
acceptable  thresholds  for  aquatic  life.  Metalimnetic  plankton  assem¬ 
blages  may  be  entrained  into  the  mixed  layer  and  initiate  plankton 
blooms . 

Breaking  internal  waves  may  increase  metalimnetic  and  hypolimnetic 
mixing.  Many  reservoirs  are  not  completely  cleared  of  timber  when 
impounded,  so  these  trees  act  as  barriers  to  wave  movement.  The  turbu¬ 
lence  generated  through  breaking  waves  can  mix  metalimnetic  constituent 


concentrations  and  increase  diffusion  across  the  epilimnetic/ 
metalimnetic  interface.  An  anoxic  metalimnion  may  contribute  dissolved 
nutrients;  reduced  species  such  as  manganese,  iron,  and  hydrogen  sul¬ 
fide;  and  organic  compounds  to  the  mixed  layer. 

Project  operation  has  the  potential  to  generate  internal  waves  of 
varying  amplitudes  and  periods.  Hydropower  generation  may  result  in 
internal  waves  with  relatively  regular  periods  while  storm  flows  can 
result  in  internal  waves  with  relatively  large  amplitudes. 

3.3,3  Turbulence,  Turbulent  Mixing,  and  Entrainment 

In  the  vertical  dimension,  the  turbulent  mixing  and  the  resulting 
constituent  distributions  are  determined  primarily  by  density  stratifi¬ 
cation.  During  periods  of  turbulent  mixing,  the  upper  well-mixed  layer 
is  usually  isotropic  despite  settling  and  vertical  variations  in  light 
intensity.  Fee  (1976)  gives  several  examples  where  in  vivo  chlorophyll 
was  isotropic  in  the  upper  mixed  layer. 

Below  the  upper  mixed  layer,  in  the  metalimnion,  the  density  gra¬ 
dient  is  usually  strong  enough  to  inhibit  turbulent  mixing.  Measure¬ 
ments  of  temperature  microstructure  in  freshwater  lakes  (Simpson  and 
Woods  1970;  Neal,  Neshyba,  and  Denner  1971)  show  that  the  density 
structure  in  the  metalimnion  is  not  smooth  as  traditionally  thought,  but 
rather  is  made  up  of  a  series  of  steps  on  very  small  length-scales 
(e.g.,  Figure  37).  These  isotropic  layers  of  intense  mixing  separated 
by  strong  gradients  can  have  significant  ecological  effects.  For 
example,  Whitney  (1938)  found  microstratif ication  on  the  scale  of 
centimetres  based  on  transparency  in  several  inland  lakes.  Whitney  also 
was  able  to  relate  these  abrupt  changes  in  transparency  to  organic 
content  and  bacterial  counts.  More  recently.  Fee  (1976)  found  narrow 
bands  of  high  chlorophyll  concentrations  in  the  metalimnion  and 
hypolimnion  of  several  small,  clear,  well-stratified  inland  lakes. 
Thornton,  Nix,  and  Bragg  (1980)  found  narrow  bands  of  high  coliform  bac¬ 
teria  concentrations  in  DeGray  Lake,  Arkansas.  Density  stratification 
and  turbulence  can  therefore  greatly  affect  vertical  constituent 
distributions . 


TEMPERATURE,  C 


Figure  37.  Example  of  temperature 
microstructure  (after  Neal,  Neshyba, 
and  Denner  1971) 


Simultaneous  measurements  of  chlorophyll  a  (a  measure  of  phyto¬ 
plankton  biomass)  and  various  physical  parameters  have  been  used  to 
determine  the  importance  of  physical  transport  on  phytoplankton 
distributions  (Abbott  et  al.  1984).  Using  a  time  series  of  chlorophyll 
and  temperature  at  a  fixed  location  and  fixed  depth  in  the  mixed  layer, 
Platt  (1972)  found  the  power  spectra  of  chlorophyll  a  followed  a  =5/3 
law  which  is  characteristic  of  three-dimensional  isotropic  turbulence. 
Using  the  same  data,  Denman  and  Platt  (1975)  found  significant  coherence 
between  temperature  and  chlorophyll  a  at  each  depth  but  not  between 
depths.  In  summarizing  these  studies  and  others,  Platt  and  Denman 
(1975b)  concluded  that  for  length  scales  between  100  m  and  5  km,  the 
observed  variations  in  temperature  and  chlorophyll  a  resulted  primarily 
from  physical  transport  mechanisms  including  internal  waves.  For  length 
scales  less  than  100  m,  the  fluctuations  were  damped  out  by  turbulent 
diffusion.  Powell  et  al.  (1975)  arrived  at  a  similar  conclusion  using 
spectra  of  current  speeds  and  chlorophyll  a  from  the  mixed  layer  of  Lake 
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Tahoe,  Nevada.  This  limit  for  turbulent  diffusion  is  also  consistent 
with  the  critical  length  scale  of  50  m  calculated  by  Platt  and  Denman 
(1975a)  for  the  smallest  patch  of  phytoplankton  that  can  maintain  itself 
against  diffusion.  Recently,  Abbott  et  al.  (1984)  found  the  formation 
and  dynamics  of  deep  chlorophyll  maxima  were  regulated  by  the  inter¬ 
action  of  three  important  processes:  turbulent  diffusion,  nutrient 
supply  rate,  and  light  availability.  Harris  (1980)  provides  an  excel¬ 
lent  discussion  of  the  important  temporal  and  spatial  scales  in  phyto¬ 
plankton  ecology  and  interaction  of  physical  processes  and  plankton 
assemblages . 

The  case  for  a  common  physical  mechanism  controlling  the  distri¬ 
bution  of  both  temperature  and  chlorophyll  a  at  selected  length  scales 
was  further  substantiated  by  Denman  (1976)  and  Fasham  and  Pugh  (1976). 

In  both  studies,  significant  coherence  was  found  between  temperature  and 
chlorophyll  a  measured  at  the  same  depth.  The  shapes  of  the  power 
spectra  of  the  two  parameters  were  also  similar.  The  gradual  steepening 
of  the  power  spectra  at  high  frequencies  was  attributed  to  internal 
waves  by  Denman  (1976).  Denman  (1976)  also  found  that  when  the 
chlorophyll  a  variance  was  high  (i.e.,  periods  of  high  biological 
activity),  the  coherence  between  chlorophyll  and  temperature  was  low. 

To  investigate  the  vertical  structure  in  more  detail,  Denman 
(1977)  used  a  series  of  vertical  profiles  of  chlorophyll  a  and  temper¬ 
ature.  The  results  suggested  that  most  of  the  chlorophyll  variation  was 
due  to  horizontal  advection  of  patches  past  the  observation  site  and  not 
to  internal  wave  motion.  Denman  (1977)  concluded  that  for  a  period 
characterized  by  little  biological  and  physical  activity,  a  single  depth 
series  of  chlorophyll  could  depict  changes  in  total  chlorophyll  only  if 
the  depth  is  located  in  an  area  of  no  appreciable  vertical  density 
gradient . 

Closely  related  to  turbulent  diffusion  is  turbulent  entrainment 
(Section  2.4.5).  Entrainment  is  precisely  what  happens  during  the 
seasonal  stratification  cycle  in  temperate  lakes.  Examples  of  the  onset 
and  breakdown  of  stratification  influencing  phytoplankton  blooms  are 
numerous  and  well  established  (Lund  1954,  Reynolds  and  Rogers  1976, 


Kaeff  and  Knoechel  1978,  etc.)*  Influences  on  the  blooms  include 
entraining  cells  living  on  the  sediment  surface  at  fall  overturn  (Lund 
1954);  entraining  previously  produced  chlorophyll  from  the  metalimnion 
and  hypolimnion  (Fee  1976);  entraining  nutrients  from  the  metalimnion 
(Stauffer  and  Lee  1974);  and  diluting  the  food  concentrations  to  a  level 
too  low  for  successful  feeding  (Fogg  1975).  Turbulent  entrainment  also 
occurs  at  shorter  time  frames  such  as  the  passage  of  synoptic  weather 
fronts.  Gusting  winds  during  storm  events  can  increase  the  depth  of  the 
mixed  layer  and  entrain  constituents  from  an  anoxic  metalimnion  and/or 
hypolimnion. 

3.3.4  Langmuir  Circulation 

A  special  case  of  turbulent  mixing  in  the  surface  waters  is 
Langmuir  circulation  (Section  3.2.4).  Observations  of  Langmuir  cells 
indicate  that  drift  velocities  are  generally  strongest  in  the  zones  of 
convergence.  The  cells  are  asymmetrical  with  downwelling  velocities 
being  larger  than  upwelling  velocities.  Downwelling  velocities  on  the 
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order  of  5  cm/sec  (4.3  x  10  m/day)  are  typical  and  are  several  orders 
of  magnitude  larger  than  typical  settling  velocities  (i.e.,  0.01  to 
1  m/day).  The  currents  in  Langmuir  circulations  are  more  than  suffi¬ 
cient  to  keep  phytoplankton  and  other  particulate  constituents  (i.e., 
bacteria)  suspended  (Denman  and  Gargett  1983). 

Langmuir  cells  can  transport  phytoplankton,  bacteria,  and  other 
constituents  into  and  out  of  the  euphotic  zone.  This  transport  can  pass 
plankton  and  other  cells  through  microstratification  layers  of  higher 
nutrient  or  organic  concentrations.  Since  nutrient  uptake  is  usually 
much  faster  than  growth,  luxury  uptake  of  nutrients  by  plankton  as  they 
pass  through  these  microstrata  can  maintain  or  promote  plankton  growth. 

3.4  Synthesis 

The  preceding  review  of  reservoir  mixing  processes  and  the 
influence  of  mixing  on  reservoir  water  quality  indicated  that  in  order 
to  achieve  the  second  objective  of  this  study  (i.e.,  develop  a 
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mathematical  algorithm  for  one-dimensional  water  quality  models  which 
realistically  represents  all  major  mixing  processes  occurring  in  CE 
reservoirs) ,  the  mixing  algorithm  should  be  capable  ol  accurately 
predicting  the: 

a.  Onset  of  stratification. 

b.  Daily  variations  in  mixed-layer  depths  and  dynamics  of  short¬ 
term  mixing  events. 

c.  Metalimnetic  gradient. 

d.  Variable  mixing  in  the  hypolimnion. 

e.  Fall  overturn. 

f ‘.  Inverse  stratification  during  winter  months. 

g.  Effects  of  project  operation. 

An  accurate  representation  of  the  onset  of  stratification  is 
required  to  accurately  predict,  for  example,  the  depletion  of  dissolved 
oxygen  and  the  onset  of  phytoplankton  blooms.  As  previously  explained 
(Section  2.2),  reservoir  stratification  starts  at  the  bottom  of  a  lake 
and  moves  upward,  not  from  the  water  surface  downward.  The  extent  of 
metalimnetic  and  hypolimnetic  DO  depletion  may  depend  on  the  relation 
between  spring  runoff  and  the  onset  of  stratification.  If  spring  runoff 
occurs  prior  to  stratification,  the  nutrients  and  organic  load  may  be 
distributed  throughout  the  water  column  or  proceed  as  a  density 
overflow.  If  the  onset  of  stratification  has  occurred  before  the  spring 
runoff,  an  interflow  may  occur  that  can  result  in  a  significant  load  to 
the  metalimnion  or  upper  part  of  the  hypolimnion.  This  load  may  exert 
an  oxygen  demand  that  results  in  the  development  of  an  anoxic 
metalimnion  or  hypolimnion.  An  accurate  representation  of  the  onset  of 
stratification  is  also  critical  for  simulating  the  water  quality  of 
lakes  that  exhibit  weak  or  intermittent  stratification. 

Daily  variations  in  mixed-layer  depths  are  required  to  accurately 
simulate  the  entrainment  of  nutrient-rich  metalimnetic  waters  into  the 
epilimnion.  This  water,  when  released  in  the  euphotic  zone,  may  result 
in  phytoplankton  blooms.  The  depth  of  the  mixed  layer,  relative  to  the 
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depth  of  the  euphotic  zone,  can  also  impact  the  timing,  magnitude,  and 
composition  of  a  phytoplankton  bloom. 

Materials  tend  to  accumulate  in  the  metalimnion  because  density 
gradients  inhibit  mixing.  It  is  necessary  to  simulate  this  accumulation 
for  entrainment  into  the  epilimnion  and  to  predict  metalimnetic  oxygen 
minima.  An  accurate  representation  of  the  metalimnetic  gradient  is  also 
required  to  compute  withdrawal  zones  and  inflow  placement. 

Variable  mixing  in  the  hypolimnion  is  required  to  accurately  simu¬ 
late  the  depletion  of  dissolved  oxygen  and  material  exchanges  across  the 
sediment/water  interface.  It  is  also  required  to  predict  hypolimnetic 
temperatures  required  to  meet  downstream  temperature  objectives. 

Although  little  is  known  about  hypolimnetic  mixing,  it  is  known  that 
mixing  levels  increase  with  hydrometeorological  forcing  and  use  of 
bottom  gates. 

With  fall  overturn,  complete  vertical  mixing  returns  and  water 
quality  problems  associated  with  stratification  are  eliminated.  If  a 
significant  oxygen  demand  builds  up  during  the  stratified  period,  oxygen 
levels  may  be  depressed  for  a  few  days  at  overturn  to  satisfy  the  exist¬ 
ing  oxygen  demand.  Fall  overturn  also  signifies  the  time  when  it  is  no 
longer  possible  to  control  release  quality  (e.g.,  temperature)  through 
project  operation  (i.e.,  selective  withdrawal). 

Many  CE  reservoirs  are  located  in  cold  climates  and  experience 
water  quality  problems  during  winter  months  when  ice  and  snow  covers 
isolate  the  lakes  from  exchanges  across  the  air/water  interface.  During 
this  period  the  lakes  can  become  inversely  stratified,  and  mixing  is 
dominated  by  inflow  and  outflow  processes. 

Project  operation  influences  both  reservoir  mixing  processes  and 
water  quality.  Bottom  withdrawal  can  promote  interflows  and  underflows, 
increase  hypolimnetic  mixing  and  mass  transfer  across  the  sediment/water 
interface,  and  decrease  the  buildup  of  anoxic  constituents  in  the  hypo¬ 
limnion.  Bottom  withdrawal  can  also  decrease  the  thermal  stability  of 
the  reservoir  by  depleting  the  colder  hypolimnetic  water  and  hastening 
fall  overturn.  Since  phytoplankton  succession  is  dictated,  in  part,  by 
the  stratification  pattern  in  the  reservoir,  the  successional  pattern 
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can  be  altered  by  project  operation.  Superimposing  peaking  hydropower 
generation  and  bottom  withdrawal  can  create  additional  mixing  through 
fluctuating  water  levels,  Kelvin-Helmholtz  instabilities,  and  epimeta- 
limnetic  and  metahypolimnetic  interfacial  shear  due  to  return  currents 
established  by  project  operation. 
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SECTION  4.  REVIEW  OF  ONE-DIMENSIONAL  PREDICTIVE  TECHNIQUES 


4.1  Introduction 


Since  the  early  to  mid  1960's,  one-dimensional  lake  and  reservoir 
models  have  been  proven  to  be  effective  tools  for  analyzing  in-lake  and 
downstream  temperature  and  water  quality  problems  because  temperature 
and  many  water  quality  parameters  tend  to  vary  more  along  a  vertical 
distance  of  tens  of  metres  than  along  a  horizontal  distance  of  thousands 
of  metres.  One-dimensional  models  solve  a  set  of  one-dimensional  con¬ 
servation  equations  for  heat  and  mass.  The  major  difficulty  in  solving 
these  equations  is  finding  expressions  for  the  turbulent  fluxes  for  heat 
and  mass. 

Following  Niiler  and  Kraus  (1977),  one-dimensional  material  trans¬ 
port  models  for  reservoirs  can  be  grouped  into  four  types  depending  on 
how  the  turbulent  fluxes  are  expressed.  These  types  include  determin¬ 
istic  solutions,  turbulence  closure  models,  eddy  coefficient  models,  and 
mixed-layer  models.  Deterministic  solutions  directly  determine  the 
fluctuating  velocities  from  the  primitive  equations  and  thereby  avoid 
the  problem  of  specifying  the  turbulent  fluxes.  Since  this  approach 
requires  a  very  fine  space  scale  and  correspondingly  high  time  reso¬ 
lution,  it  is  too  expensive  and  time  consuming  for  practical  problems. 

An  example  of  its  use  is  found  in  Deardorff  (1970).  Turbulence  closure 
models  express  the  turbulent  fluxes  in  terms  of  higher  order  moments 
(i.e.,  averaged  triple  products  of  the  fluctuating  quantities)  which  in 
turn  must  be  parameterized  in  terms  of  empirical  coefficients  and  com¬ 
putable  (bulk)  quantities  or  be  represented  by  another  set  of  equations 
involving  fourth-order  moments,  Mellor  and  Yamada  (1974)  describe  this 
process  as  a  hierarchy  of  turbulent  closure  models  with  the  classical 
eddy  coefficient  models  at  the  lowest  level.  With  the  exception  of  this 
lowest  level,  the  equations  solved  in  these  models  are  cumbersome  and 
difficult  to  solve.  Eddy  coefficient  models  assume  the  turbulent  fluxes 
can  be  expressed  by  the  gradient  of  the  transported  quantity  multiplied 
by  an  eddy  diffusion  coefficient  which  is  a  complicated  function  of 
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space  and  local  stability.  Mixed-layer  models  assume,  as  implied  by  the 
name,  that  the  upper  layer  is  well  mixed.  This  assumption  permits  the 
vertical  integration  and  solution  of  the  turbulent  energy  equation  in 
terms  of  the  upper  mixed-layer  depth.  According  to  Niiler  and  Kraus 
(1977),  this  approach  is  probably  the  most  effective  tool  for  prediction 
of  upper  ocean  temperatures  and  heat  storage. 

In  this  section,  the  one-dimensional  assumption  and  its  implica¬ 
tions  are  examined,  and  existing  eddy  coefficient  and  mixed-layer  models 
are  reviewed.  Deterministic  solutions  and  turbulence  closure  models 
were  not  considered  cost-effective  tools  for  the  solution  of  practical 
problems . 


4.2  One-Dimensional  Assumption 

Vertical  one-dimensiouality  assumes  that  forcing  across  the  air/ 
water  interface  and  vertical  variations  (gradients)  play  a  dominant  role 
with  transverse  and  longitudinal  variations  playing  a  secondary  role. 

The  one-dimensional  assumption  therefore  assumes  that  reservoirs  are 
represented  by  a  vertical  series  of  horizontal  slices  (Figure  38).  Each 
horizontal  slice  is  assumed  to  be  well  mixed  laterally  and  longitudi¬ 
nally.  Inputs  to  a  horizontal  layer  are  assumed  to  be  instantaneously 
mixed  throughout  the  layer.  Since  vertical  density  gradients  (i.e., 
strong  stratification)  inhibit  vertical  motions  (i.e.,  act  to  keep  the 
isotherms  horizontal) ,  the  use  of  a  one-dimensional  assumption  for 
stratified  systems  is  often  justified. 

Ford  and  Thornton  (1979)  analyzed  the  time  and  length  scales  char¬ 
acterizing  the  hydrodynamics,  chemistry,  and  biology  of  lakes  and  showed 
that  there  is  both  an  upper  and  lower  bound  to  the  lake  size  that  can  be 
characterized  by  a  one-dimensional  model.  The  lower  bound  was  set  at 
horizontal  dimensions  on  the  order  of  0.05  to  0.5  km  to  minimize  hori¬ 
zontal  gradients  caused  by  near-field  effects  such  as  local  runoff  and 
upwelling.  The  upper  bound  was  set  at  horizontal  dimensions  on  the 
order  of  10  to  1.0  km  to  minimize  horizontal  variations  caused  by  vari¬ 
ations  in  the  synoptic  forcing  functions. 
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TRIBUTARY 

INFLOW 


Figure  38.  Vertical  one-dimensional 
representation 

When  evaluating  the  appropriateness  of  the  one-dimensional  assump¬ 
tion  for  reservoirs,  two  physical  factors  must  be  specifically  consid¬ 
ered  in  addition  to  the  ideas  presented  by  Ford  and  Thornton  (1979). 
These  factors  are  the  wind  setup  and  the  horizontal  penetration  of  the 
inflow  (i.e.,  plunge  point).  When  the  wind  blows  across  the  water  sur¬ 
face,  the  resulting  shear  stress  causes  the  isotherms  to  tilt 
(Section  3.2.2).  If  the  wind  is  sufficiently  strong,  metalimnetic  and 
hypolimnetic  waters  may  become  exposed  to  the  surface  (upwelling)  and 
become  mixed  with  epilimnetic  waters.  The  magnitude  of  the  wind  setup 
at  the  water  surface  can  be  estimated  by  comparing  the  applied  wind 
shear  force 
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t  A  =  p  C  W^A  =  p  w^A 
ss  ad  s  w  *  s 


where 
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x  =  shear  stress,  kg/(m-sec  ) 

S  2 
A  *  surface  area,  m 

S  3 

p  -  density  of  air,  1.177  kg/m 

3. 

W  =  wind  speed,  m/sec 
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p  =  density  of  water,  kg/m 
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with  the  pressure  force  (potential  energy)  that  builds  up 
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This  balance  assumes  no  bottom  shear.  The  slope  of  the  water  surface  is 
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and  the  slope  of  the  thermocline 
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If  (dz/dx)  is  approximated  by  AD/L,  Imberger's  Wedderburn  number  (We)  is 
obtained 
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The  Wedderburn  number  therefore  compares  the  slope  of  the  interface  with 
wind  forcing.  It  can  be  used  to  determine  if  upwelling  is  significant. 
If  We  »  1,  upwelling  is  not  sufficient  and  the  one-dimensional  assump¬ 
tion  is  valid. 

The  plunge  point  location  provides  an  indication  of  where  the 
riverine  zone  stops  and  the  lake  begins  (Section  3.2.5).  If  the  river¬ 
ine  zone  extends  a  significant  distance  into  the  reservoir,  the  reser¬ 
voir  is  dominated  by  advective  forces  and  cannot  be  considered  to  be 
vertically  one-dimensional.  Following  the  recommendations  of  Ford  and 
Johnson  (1983),  the  depth  of  the  plunge  point  can  be  determined  by 
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where 


h 

P 


B 

c 


hydraulic  plunge  depth,  m 

3 

inflow  rate,  m  / sec 

width  of  zone  of  conveyance,  m 

3 

density  difference  between  inflow  and  reservoir,  kg/m 


If  h  /Z  ,  where  Z  is  the  maximum  reservoir  depth,  is  less  than  1/3, 
p  m  m 

the  reservoir  is  one-dimensional. 


4.3  Eddy  Coefficient  (Diffusion)  Models 


4.3.1  Description 

Eddy  diffusion  models  were  initially  used  in  the  mid-1960’s  to 
predict  temperature  changes  in  deep,  stratified  reservoirs.  They  have 
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been  routinely  used  by  the  CE  since  the  late  1960's  and  form  the  basis 

for  many  one-dimensional  water  quality/ecological  models,  Eddy 
diffusion  models  are  based  on  the  vertical  one-dimensional  thermal 
energy  equation; 


f +  iky  fe  <v>  ■  iky  fc  <K» +  V  A(z)  f 


+  T7£)  "  (z)T  -  u  (z)t] 

A(*i  <  J  J  0  J 


pwCpA(z)  3z  [Hz  A(z)] 


where 


UjU) 

T. 

Uo(Z) 


=  water  temperature,  °C 
=  time,  sec 

r* 

=  horizontal  area  of  reservoir  at  elevation  z,  m 
■  vertical  elevation,  m 

3 

=  vertical  flow  rate,  m  /sec. 

2 

=  molecular  diffusion  coefficient,  m  /sec 

2 

=  global  vertical  diffusion  coefficient,  m  /sec 

*=  reservoir  width  at  elevation  z,  m 

“  inflow  velocity  distribution,  m/sec 

=  inflow  temperature,  °C 

“  outflow  velocity  distribution,  m/sec 

3 

density  of  water,  kg/m 

=  heat  capacity  of  water,  J/(kg-°C) 

2 

r  heat  flux  at  elevation  z,  W/m 


Term  (1)  on  the  left  side  is  the  time  rate  of  change  of  tempera¬ 
ture.  Term  (2)  is  the  vertical  advective  transport  term.  It  represents 
the  flow  between  internal  elements  that  is  required  to  balance  inflows 
and  outflows  and  maintain  continuity.  If  a  Lagrangian  or  variable  layer 


scheme  is  used  (Section  5.4.3),  there  is  no  flow  between  elements 
(Qv  =  0),  and  term  (2)  is  not  required.  Term  (3)  (the  first  term  on 
right  side  of  Equation  29)  is  the  diffusion  term.  The  molecular  diffu¬ 
sion  coefficient  Km  ,  is  usually  neglected  or  incorporated  into  the 

turbulent  diffusion  coefficient  K  .  Various  formulations  for  K  are 

z  z 

discussed  in  Section  4.3.2.  Term  (4)  of  Equation  29  represents  the  heat 
input  from  tributaries  and  heat  loss  from  outflows.  The  inflow  and  out¬ 
flow  velocity  distributions  are  usually  based  on  empirical  relationships 
developed  from  laboratory  experiments  (e.g. ,  Bohan  and  Grace  1973). 

Many  diffusion  models  developed  for  natural  lakes  (e.g.,  Sundaram  and 
Rehm  1973;  Henderson-Sellers  1976;  Walters,  Carey,  and  Winter  .1978) 
ignore  the  fourth  term.  For  many  CE  reservoirs  characterized  by  short 
residence  times  (i.e.,  less  than  60  days),  this  term  dominates 
Equation  29.  Term  (5),  the  last  term  in  Equation  29,  represents  the 
local  heat  source  which  includes  the  internal  heating  from  solar 
radiation.  The  various  components  of  this  term  were  discussed  in 
Section  3.2.1. 

Equation  29  can  be  solved  by  various  computational  techniques 
knowing  the  initial  conditions,  boundary  conditions  (i.e.,  heat  transfer 
at  air/water  interface,  inflow  and  outflow  velocity  distributions,  heat 
transfer  at  sediment/water  interface),  and  reservoir  geometry  (i.e., 

A(z)  and  B(z)).  In  addition  to  these  requirements,  the  only  major  dif¬ 
ficulties  remaining  to  solve  Equation  29  are  the  specifications  of  the 
internal  absorption  of  solar  radiation  and  the  eddy  diffusion 

coefficient  Kz  .  As  indicated  in  Section  2.2.2,  the  specification  of 
requires  knowledge  of  both  dissolved  and  particulate  materials  in 
the  water.  Although  thermal  simulations  in  some  systems  are  sensitive 
to  the  attenuation  of  shortwave  radiation,  it  will  not  be  discussed 
further.  Reviews  on  light  penetration  can  be  found  in  Williams  et  al. 
(1981) . 

4.3.2  Eddy  Coefficient  Formulations 

Field  measurements.  Since  the  eddy  diffusion  coefficient  includes 
the  dynamics  of  all  mixing  processes,  it  is  a  complicated  function  of 

8C 


m  "  # 


•  .  m  .. 


time,  space,  and  local  stability.  It  is  therefore  informative  to  review 
results  from  field  measurements  of  eddy  diffusion  coefficients  prior  to 
discussing  the  various  formulations  for  .  According  to  Imberger  and 

Hamblin  (1982),  there  are  no  direct  measurements  of  local  values  of 
vertical  diffusion  coefficients  in  either  lakes  or  oceans,  only  basin 
averaged  values  derived  from  global  thermal  budgets  or  tracer 

budgets.  Since  this  study  is  concerned  only  with  one-dimensional 
formulations,  global-  or  basin-averaged  coefficients  are  consistent  with 
study  assumptions. 

In  the  global  thermal  budget  approach,  estimates  of  K  are  back- 

z 

calculated  from  Equation  29,  temperature  data,  and  information  on  the 

internal  absorption  of  solar  radiation  Hz  ,  assuming  no  vertical  and 

horizontal  advective  transport  (i.e.,  terms  (2)  and  (4)  in  Equation  29 

are  assumed  to  be  0) .  This  assumption  is  usually  valid  for  lakes  and 

reservoirs  with  long  residence  times.  If  Equation  29  is  integrated  down 

the  water  column  from  depth  z  to  the  maximum  depth  Z  ,  then 

m 

“■•/■( B  ".)/(-«) 

and  K?  can  be  determined  either  graphically  or  with  the  aid  of  a  com¬ 
puter.  Typical  results  are  shown  in  Figure  39,  which  illustrates 
decreased  values  for  Kz  in  the  metalimnion  or  regions  of  temperature 
gradient.  The  values  computed  by  the  method  represent  an  average 

value  over  the  period  between  temperature  profiles.  Figure  40  illus¬ 
trates  that  different  values  can  be  obtained  if  different  integration 
periods  are  used.  The  smaller  the  period,  the  more  variation  in  K 

z 

Computed  Kz  values  may  vary  over  several  orders  of  magnitude  from  one 
day  to  the  next.  As  discussed  in  Jassby  and  Powell  (1975),  it  is  essen¬ 
tial  that  the  internal  absorption  of  solar  radiation  be  considered  when 

back-calculating  K  .  For  this  reason  and  others,  K  should  not  be 

z  z 

computed  for  the  epilimnion  using  this  method. 
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Figure  39.  Vertical  variation  in  eddy 
diffusion  coefficients 


Eddy  diffusion  coefficients  can  also  be  computed  from  passive 
tracer  budgets  using  the  one-dimensional,  unsteady  diffusion  equation 
(i.e.,  a  modified  version  of  Equation  29)  and/or  the  method  of  moments 
(Fischer  et  al.  1979).  Examples  of  tracer  studies  include  Kullenburg, 
Murthy,  and  Westerberg  (1973)  using  dye;  Imboden  et  al.  (1977)  using 
tritium;  Imboden  and  Emerson  (1978)  using  radon  and  phosphorus,  and  Quay 
et  al.  (1980)  using  tritium.  The  K values  computed  from  these  stud¬ 
ies  represent  the  time  scale  over  which  the  study  took  place.  In  many 
instances  this  period  is  shorter  than  those  computed  by  the  thermal 
budget  approach. 
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Figure  40.  Variation  in  computed  eddy 
diffusion  coefficients  with  integration 
time,  McCarrons  Lake,  Minnesota 


Field  studies  using  the  thermal  budget  approach,  tracers,  or  a 

combination  of  both  substantiate  the  following  conclusions  concerning 

vertical  diffusion  coefficients.  First,  vertical  dif fusivities  of  heat 

-7  2 

range  from  molecular  (i.e.,  1.42  *  10  m  /sec)  up  to  values  of 
-4  2 

10  m  /sec  (Figure  41).  Th< y  are  several  orders  of  magnitude  less  than 
horizontal  eddy  diffusion  coefficients  but  greater  than  molecular  diffu¬ 
sion  coefficients  for  gases  and  salts. 

Second,  vertical  diffusion  coefficients  for  heat  and  mass  are 
about  the  same  in  the  hypoliranion  but  differ  in  the  metalimnion.  Quay 
et  al.  (1980)  found  that  heat  diffuses  faster  than  mass  in  the  metalim¬ 
nion,  indicating  molecular  diffusion  is  important  in  determining  the 
rate  of  vertical  transport  in  the  metalimnion.  As  shown  in  Table  2,  the 
molecular  diffusion  coefficient  for  heat  and  mass  differ  significantly. 

Third,  the  diffusivity  is  usually  higher  during  periods  of  strong 
winds  and  large  inflows  or  outflows.  Bengtsson  (1978)  found  an  almost 
linear  increase  in  log  with  wind  speed  (Figure  42). 

Fourth,  the  diffusivity  decreases  with  increasing  stability  (i.e., 
stratification)  according  to  the  following  relationship 
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Figure  41.  Comparative  ranges  for  diffusion 
coefficients  (after  Lerman  1971) 
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Figure  42.  Increase  in  K  with  wind 
speed  (after  Bengtsson  1978) 


K  =  a(N2)_n 
2 


(31) 


where 

a  *  proportionality  coefficient,  units  vary 
N2  =  ^  ~  =  buoyancy  frequency,  sec  2 

P  o  Z 

n  =  coefficient,  dimensionless 

Values  for  n  vary  from  0.25  (Hutchinson  1941)  to  almost  2.0  (Blanton 

1973).  Using  dimensional  analysis,  Welander  (1968)  attriDuted  different 

coefficients  to  different  physical  processes.  For  example,  n  =  0.5  is 

a  local  shear  process  and  n  =  1  is  for  a  cascade  process.  Figure  43 

2 

illustrates  the  variation  in  K  with  N 

z 

Fifth,  diffusion  coefficients  generally  increase  with  lake  size. 
Ward  (1977)  related  K to  the  surface  area  of  lakes 

-9  1  /  2 

K  =  3.3  x  10  A  '  (32) 

z  s 
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Figure  43.  Variation  in  K  with  stratification 
2  2 

stability  N  (after  Quay  at  al.  1980) 


where  A  =  surface  area,  m  . 
s 


This  relationship  is  consistent  with  the  data  published  by  Ford 
(1978)  (Figure  44). 


Figure  44.  Variation  in  K  with 
lake  surface  area 

Formulations.  The  specification  of  K  varies  from  model  to 
model.  If  turbulent  diffusion  Is  assumed  to  be  much  larger  than 
molecular  diffusion,  molecular  diffusion  can  be  incorporated  into  the 
turbulent  term  without  any  loss  In  precision.  Some  early  thermal  models 
(e.g.,  MIT  model  (Huber  and  Harleman  1968);  Eiker  Model  (USAE  District, 
Baltimore  1974))  considered  the  turbulent  diffusivity  to  be  constant 
with  depth;  others  considered  it  to  be  dependent  on  the  density  struc¬ 
ture  and,  therefore,  depth  dependent. 

Historically,  the  steady-state  model  of  Munk  and  Anderson  (1948) 
was  probably  the  first  to  provide  insight  into  the  interplay  of  turbu¬ 
lence,  stratification,  and  heat  flow.  They  assumed  the  eddy  diffusion 
coefficient  was  related  to  the  Richardson  number  by 
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(33) 


K  =  K  (1  +  a  Ri) 
z  zo 


where 

2 

K  =  eddy  diffusion  coefficient  at  neutral  stability,  m  / sec 
a  =  dimensionless  coefficient 
Ri  =  Richardson  number,  dimensionless 


Since  then,  variations  of  this  formulation  have  been  used  by  Sundaram 

and  Rehm  (1971,  1973);  Henderson-Sellers  (1976);  Walters,  Carey,  and 

Winter  (1978);  and  others.  Specific  differences  between  formulations 

centered  around  the  specification  of  K  ,  Ri  ,  and  the  power  (i.e., 

zo 

-3/2).  Some  specific  models  (e.g.,  Sundaram  and  Rehm)  modify  the  form¬ 
ulation  after  stratification  forms  to  account  for  the  smaller  eddy 
dif fusivities  that  are  normally  found  in  the  hypolimnion  of  stratified 
lakes . 

Henderson-Sellers  (1976)  reviewed  five  expressions  for  K  : 

zo 


K 


zo 


Pr 


wJ 


3U 

3z 
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zo 


p? k  u* 


(z„  -  z) 


zo 


cv. 
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1,2,  ,2  3U 

P?  k  <ks  -  te 


K 

zo 


_L  k2 

Pr  K 


3U 

3z 


32U 


3z 


(34) 

(35) 

(36) 

(37) 

(38) 


where 

Pr  =  Prandtl  number  =  1  for  water,  dimensionless 
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w*  =  shear  velocity,  m/sec 
8  U 

*  vertical  velocity  gradient,  (m/sec)/m 

k  =  von  Karman's  constant  ('■0.4) 

z  =  elevation  of  water  surface,  m 
s 

z  =  elevation,  m 
c  =  empirical  coefficient,  m 

and  recommended  Equation  34.  The  major  difficulty  with  Equations  34, 

37,  and  38  is  that  they  require  a  priori  knowledge  of  the  current  struc¬ 
ture.  Although  it  is  possible  to  assume  a  simple  characteristic  veloc¬ 
ity  profile  such  as  parabolic  or  exponential  under  isothermal,  steady- 
state  conditions,  stratification  significantly  modifies  these  profiles. 
In  addition,  the  impact  on  the  current  profile  of  the  time-varying 
inflow  density  currents,  withdrawal  zones-  wind  vectors,  and  complicated 
reservoir  geometry  is  totallj'  unknown.  The  practicality  of  formulations 
for  K  which  require  the  velocity  gradient  is  therefore  questionable. 
Physically,  Equation  35  is  questionable,  since  increases  with 

depth,  which  is  contrary  to  field  observations.  Equation  36  is 
desirable  since  the  Kzq  increases  with  wind  shear  velocity  in  water 

but  undesirable  since  K  is  constant  with  depth  and  c  is  not 

zo  r 

dimensionless . 

The  form  of  the  Richardson  number  (Ri)  in  Equation  33  can  be  taken 
as 


or 


(39) 


(40) 
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where 


g  «  acceleration  due  to  gravity,  m/sec 
~  =  local  density  gradient 

As  with  the  formulations  for  K  ,  Equation  40  is  preferable  over 

zo 

Equation  39  since  the  velocity  gradient  is  not  known. 

Substitution  of  Equation  36  and  40  into  Equation  33  yields 


K  =  - 
z 

_  £  /l£\ 

p  Uz/ 

1b 

< 

i  t*  a 

2 

w* 

► 

L  <zs  -  2)2  J 

(41) 


If  a  lake  does  not  stratify  (i.e.,  Ri  =  0),  Equation  41  reduces  to 

K  =  cw.  (42) 

z  * 

and  if  a  lake  stratifies  such  that  Ri  »  1,  Equation  41  reduces  to 


K 


where  c* 


in  met  t  is. 


(43) 


Ozmidov  (1965)  too*  a  different  approach  to  formulating  an  edd) 
diffusion  coefficient  .  lie  assumed  that  mixing  was  done  at  the  smaller 
length  scales  and  thud:  a  lxna.l  steady  mte  e\  sted  such  that  the  rate 
of  turbulent  unergy  ztnptn  a:  la  ger  length  sc,*  es  is  in  balance  with 
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rate  of  dissipation  at  smaller  scales.  Ozmidov  (1965)  used  an  equation 
similar  to  Equation  43  with  b  =  1 


where 


c  e 

Kz  =  ' — 2 
N 


=  global  vertical  diffusion  coefficient,  m  /sec 
=  dimensionless  calibration  coefficient 

2  3 

=  dissipation  rate  of  TKE,  m  /sec 
-  p~  =  bu°yancy  frequency,  sec  ^ 


If  all  of  the  TKE  is  derived  from  the  wind,  the  local  dissipation  per 
unit  volume  is  proportional  to  the  TKE  input  per  unit  surface  area 

3 

(i.e.,  a  w^)  (see  Equations  53  and  54)  divided  by  the  water  column  depth 
z  .  That  is. 


and  Equation  45  is  similar  to  the  numerator  of  Equation  43  if  b  -  I  . 

Other  formulations  for  eddy  diffusion  coefficients  have  been 
proposed,  but  they  are  all  variations  of  the  above  relationships  and 
therefore  will  not  be  discussed. 


4,3.3  Assumptions  and  Limitations 

The  major  assumption  of  eddy  diffusion  models  is  that  the  effects 
of  all  mixing  processes  can  be  combined  into  a  single  diffusion  coeffi¬ 
cient  K?  ,  which  is  coupled  to  the  mean  concentration  gradient  9C/9z 
(see  Section  2.4.4).  If  the  concentration  gradient  becomes  small  or 
goes  to  zero  (i.e.,  as  in  the  epilimnion) ,  the  diffusion  coefficient 
must  be  made  infinitely  large  to  compensate  for  the  decrease  in  gradient 
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and  still  result  in  a  finite  value  of  the  flux,  K  3C/3z.  The  method, 

z 

therefore,  breaks  down  when  vertical  transport  is  under  way.  According 
to  Tennekes  and  Lumley  (1972),  the  concept  makes  sense  only  for  flows 
with  a  single  characteristic  scale.  Since  mixing  in  reservoirs  is  char¬ 
acterized  by  many  scales  of  motion,  the  validity  of  the  concept  is 
questionable.  Since  turbulence  may  be  decoupled  from  the  mean  gradient, 
the  concept  is  physically  meaningless  (Zeman  1981). 

From  a  practical  viewpoint,  the  one  major  limitation  of  the  vari¬ 
ous  empirical  formulations  for  the  K^  is  lack  of  knowledge  concerning 
the  internal  current  structure  (i.e.,  3U/3z)  in  reservoirs.  Since  the 
one-dimensional  models  under  consideration  in  this  study  do  not  and  are 
not  capable  of  simulating  the  internal  current  structure  of  reservoirs 
(i.e.,  that  is  not  their  purpose),  formulations  for  Kz  that  include 
3U/3z  are  not  practically  viable  alternatives. 


4.4  Mixed-Layer  Models 

4.4.1  Background 

In  reservoirs,  mixing  is  primarily  caused  by  turbulence  (Sections 
2.4.4  and  3.2.3).  Since  turbulence  is  highly  dissipative,  it  requires  a 
constant  source  of  energy  to  be  maintained.  In  addition,  energy  is  also 
expended  when  mixing  occurs  (Section  2.3).  It  is  therefore  logical  that 
a  turbulent  kinetic  energy  (TKE)  budget  be  used  to  analyze  mixing.  TKE 
budgets  have  been  used  by  atmospheric  scientists  to  study  entrainment  at 
the  inversion  base  (e.g.,  Tennekes  1973;  Stull  1976a, b;  Zeman  and 
Tennekes  .1977;  Deardorff  1979;  Mahrt  1979),  by  oceanographers  to 
investigate  entrainment  into  the  upper  mixed  layer  (e.g.,  Kraus  and 
Turner  1967;  Turner  and  Kraus  1967;  Pollard,  Rhines,  and  Thompson  1973; 
Niiler  1975;  Niiler  and  Kraus  1977;  Garnich  and  Kitaigorodskll  1977, 
1978),  and  by  physical  limnologists  to  study  entrainment  in  lakes  and 
reservoirs  (e.g.,  Stefan  and  Ford  1975a, b;  Hurley-Octavio,  Jirka,  and 
Harleman  1977;  Imberger  et  al.  1978;  Bloss  and  Harleman  1980;  Ford  and 
Stefan  1980b;  Imberger  and  Patterson  1981).  Atmospheric  scientists 
generally  perform  the  energy  budget  at  the  density  interface,  while 


oceanographers  and  engineers  generally  employ  a  vertically  integrated 
energy  budget.  These  models  are  therefore  sometimes  referred  to  as 
integral  energy  models  or  mixed-layer  models.  According  to  Tennekes  and 
Driedonks  (1980),  there  are  no  substantial  differences  in  these  two 
approaches.  The  differences  deal  primarily  with  notation. 

Mixed-layer  models  have  their  origin  with  the  laboratory  experi¬ 
ments  of  Rouse  and  Dodu  (1955).  They  studied  mixing  across  a  density 
interface  by  generating  turbulence  in  the  upper  layer  of  a  two-layered 
density-stratified  fluid  using  an  oscillating  grid  (Figure  45).  The 
turbulence  formed  a  well-mixed  layer,  bounded  by  a  sharp  interface  which 
moved  away  from  the  stirrer  as  fluid  was  entrained  from  the  underlying 
quiescent  layer  into  the  upper  well-mixed  layer.  Since  then,  other 
experiments  employing  oscillating  grids  (e.g.,  Cromwell  1960;  Bouvard 
and  Dumas  1967;  Turner  1968;  Brush  1970;  Wolanski  1972;  Crapper  1973; 
Crapper  and  Linden  1974;  Linden  1973,  1975;  Thompson  and  Turner  1975; 
Wolanski  and  Brush  1975;  Hopfinger  and  Toly  1976),  moving  screens  at  the 
water  surface  (e.g.,  Kato  and  Phillips  1969;  Kantha,  Phillips,  end  Azad 
1977),  wind  (e.g,,  Wu  1973),  and  oppositely  directed  jets  (e.g.,  Moore 
and  Long  1971),  have  indicated: 

a.  The  turbulence  produces  and  maintains  a  well-mixed  layer, 

b.  The  turbulence  sharpens  the  interface. 

£.  The  turbulence  causes  mixing  across  the  interface. 

d.  The  entrainment  velocity  or  advancement  speed  of  the  interface 
is  a  function  of  the  Richardson  number. 

e.  Mixing  takes  place  largely  through  a  process  which  looks  like 
the  intermittent  breaking  of  steep-faced  internal  waves  which 
tend  to  thicken  the  interface,  followed  by  the  sweeping  away 
of  this  fluid  by  the  stirring  in  the  layers,  which  sharpens 
the  interface  again. 

f.  .  The  turbulence  reflects  the  characteristics  of  the  generation 

mechanism  (e.g.,  grid),  and  therefore  the  scale  of  the 
generated  turbulence  is  smaller  than  the  mixed-layer  depth. 

£.  The  TKE  generated  by  stirring  mechanisms  decays  with  distance 
from  the  mechanism. 

The  results  from  these  experiments  are  summarized  in  Turner  (1973); 
Sherman,  Imberger,  and  Corcos  (1978);  and  Pederson  (1980). 
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Figure  45.  Schematic 
of  an  oscillating  grid 
laboratory  experiment 


4.2  Description 

Results  from  the  experiments  described  in  Section  4.4.1  and  others 
are  used  to  parameterize  and  solve  the  TKE  equation,  which  is  the  basis 
of  mixed-layer  models.  The  TKE  equation  is  derived  by  scalar  multipli¬ 
cation  of  the  momentum  equations  with  the  turbulent  velocity  field  and 
assuming;  (a)  the  Boussinesq  approximation;  (b)  the  main  flow  is 
horizontal  and  in  the  direction  of  U  ;  and  (c)  the  turbulence  is  homo¬ 
geneous  in  the  horizontal  plane  (see,  e.g.,  Van  Mieghem  1973,  Ford  1976, 
and  Phillips  1977): 


q2  +  --  (u  (  ^—  +  q2))  =  -u  u  t-~  -  u  p ’  -  e  (46) 

at  M  az  pQ  x  z  9z  p0  z 

(1)  (2)  (3)  (4)  (5) 

where 

q2  -  (u2  +  u2  +  u2)/2  =  TKE  per  unit  mass 
x  y  z 

u  ,  u  ,  u  *  turbulent  fluctuations  of  the  horizontal  (x,  y)  and 
A  y  x  vertical  (z)  velocity  components,  m/sec 
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p*  ■  turbulent  fluctuations  of  pressure.  Pa 

2 

p’  D  turbulent  fluctuations  of  density,  kg/m 

2 

p  “  mean  density,  kg/m 

°  2  3 

e  ■  rate  of  dissipation  of  TKE,  m  /sec 

■  mean  horizontal  velocity,  m/sec 

t  ■  time,  sec 

z  =  vertical  coordinate,  m 

In  Equation  46,  term  (1)  represents  the  temporal  change  in  TKE;  term  (2) 
represents  the  redistribution  in  space  of  the  TKE  by  turbulence; 
term  (3)  represents  the  rate  of  transfer  of  TKE  from  the  mean  flow  by 
the  Reynolds  stresses;  term  (4)  represents  the  gain  or  loss  of  TKE  due 
to  the  release  or  increase  of  potential  energy  of  the  mean  density  (tem¬ 
perature)  field;  and  term  (5)  represents  the  rate  of  dissipation  of  TKE. 

In  its  present  form.  Equation  46  cannot  be  solved  directly  for  any 
parameter  of  practical  Interest.  The  equation  must  therefore  be  parame¬ 
terized  in  terms  of  variables  of  interest  and  practical  significance. 
Following  the  parameterization  scheme  of  Tennekes  and  Driedonks  (1980), 
Equation  46  becomes 
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where 

Ct’  Cf’  Cd  =  constants  determined  from  experimental  data 
o  =  velocity  scale,  m/sec 
h  *  depth  of  mixed  layer,  m 

w  =  =  entrainment  velocity,  m/sec 

6  dt 

Early  mixed-layer  models  (i.e.,  Kraus  and  Turner  1967;  Stefan  and 
Ford  1975a,b;  Hurley-Octavio ,  Jirka,  and  Harleman  1977;  Imberger  et  al. 
1970)  assumed  the  temporal  and  shear  terms  (terms  (1)  and  (3)  in  Equa¬ 
tion  46)  were  zero  and  that  the  local  dissipation  equaled  local 


production  such  that  Equation  47  reduced  to: 
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when  o  was  assumed  to  scale  with  .  This  dependence  of  the 

entrainment  speed  on  the  inverse  of  the  Richardson  number  was  found  in 
many  of  the  laboratory  experiments  previously  discussed  in  this  section. 
Models  based  on  Equation  48  did  an  excellent  job  of  simulating  mixed- 
layer  dynamics  during  periods  of  strong  stratification  but  overpredicted 
mixed-layer  depth  during  periods  of  weak  stratification  (Bloss  and 
Harleman  1980) . 

Bloss  and  Harleman  (1980)  corrected  this  deficiency  by  using  the 
parameterization  of  Zilitinkevich  (1979)  for  the  transient  term 
(term  (1)  in  Equation  46).  Assuming  shear  was  still  small  and  using  the 
constant  from  Zeman  and  Tennekes  (1977),  Bloss  and  Harleman  (1980) 
expressed  Equation  47  as  a  function  of  Richardson  number. 
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This  function  (Figure  46)  is  actually  an  efficiency  factor  for  the  con¬ 
version  of  TKE  into  potential  energy.  This  modification  significantly 
improved  model  predictions  during  periods  of  weak  stratification. 

Another  deficiency  of  mixed-layer  models  was  noted  by  Imberger  and 
Patterson  (1981).  They  observed  that  mixed-layer  models  always  cause 
the  metalimnetic  gradient  to  sharpen,  while  field  observations  indicated 
times  when  the  gradient  was  more  diffuse.  This  prompted  Imberger  and 
Patterson  to  include  shear  (term  (3)  in  Equation  46)  in  their 
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Figure  46.  Efficiency  factor  for  mixing 
(after  Bloss  and  Harleman  1980) 

formulation.  This  improvement  requires  computing  the  slope  of 
thermocline  under  wind  forcing. 

4.4.3  Assumptions  and  Limitations 

TKE  or  mixed-layer  models  assume  that  a  well-mixed  layer  overlies 
a  region  of  density  gradient.  The  depth  of  the  upper  mixed  layer  is 
determined  by  the  TKE  input  at  the  air/water  interface  and  the  retarding 
buoyancy  force  of  the  underlying,  density-srratif ied  layer.  The  models 
predict  averaged  parameters  for  the  mixed  layer  but  do  not  consider  the 
dynamics  of  the  underlying  density-stratified  layer.  Its  major  limita¬ 
tion  is,  therefore,  that  it  does  not  simulate  the  lower  density- 
stratified  layer. 


4 . 5  Summary 

The  review  of  diffusion  and  TKE  models  indicates  that  both 
approaches  have  significant  advantages  and  limitations.  More  important, 
however,  is  the  fact  that  the  two  approaches  model  different  physical 
mixing  processes  (i.e.,  diffusion  and  entrainment)  and,  since  both  of 
these  processes  are  important  reservoir  mixing  processes,  both  should  be 
included  in  the  recommended  algorithm. 
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SECTION  5.  ALGORITHM  DEVELOPMENT 


5.1  Introduction 


The  second  objective  of  this  study  was  to  develop  a  mathematical 
algorithm  that  realistically  represents  all  major  mixing  processes 
occurring  in  CE  reservoirs.  Since  this  algorithm  is  intended  to  be  used 
in  a  generalized  one-dimensional  water  quality  model  (e.g.,  CE-QUAL-R1 
(Environmental  Laboratory  1982))  to  predict  changes  In  reservoir  water 
quality  resulting  from  changes  in  hydrometeorological  conditions  and 
project  operation,  several  requirements  had  to  be  considered  during 
algorithm  development. 

First,  the  algorithm  had  to  be  generalized  with  respect  to  CE  res¬ 
ervoirs.  Sizes  of  CE  reservoirs  range  over  four  orders  of  magnitude 
(Ford  1978).  Some  reservoirs  are  small  and  round  (e.g.,  Eau  Galle  Lake, 
Wisconsin,  and  C.  J.  Brown  Lake,  Ohio)  while  others  are  large  and  den¬ 
dritic  in  shape  with  many  islands  (e.g.,  DeGray  Lake,  Arkansas).  Ele¬ 
vations  and  proximity  to  major  cities  (i.e.,  a  source  of  meteorological 
data)  also  vary.  Since  there  is  no  typical  CE  reservoir  or  location, 
the  algorithm  must  therefore  not  be  constrained  by  extensive  morpho¬ 
metric  and  hydrometeorological  data  requirements.  In  addition.  It  Is 
highly  unlikely  that  similar  mixing  processes  dominate  in  all 
reservoirs. 

Second,  the  algorithm  is  to  be  used  in  a  one-dimensional  water 
quality  model  (CE-QUAL-R1) .  The  one-dimensional  assumption  (Sec¬ 
tion  4.2)  not  only  limits  model  applicability  to  certain  types  of  lakes 
and  problems  but  also  limits  the  types  of  formulations  that  can  be  used. 
For  example,  one-dimensional  models  do  not  generally  compute  the  inter¬ 
nal  current  structure.  Water  quality  is  defined  to  include  the  kinds 
and  amounts  of  dissolved  and  suspended  matter  in  water,  the  physical- 
chemical  characteristics  of  the  water,  and  the  ecological  relations  with 
and  among  aquatic  organisms;  therefore,  the  algorithm  must  be  capable  of 
accurately  describing  the  vertical  transport  of  all  water  quality  param¬ 
eters  and  yet  not  be  sensitive  to  concentration  gradients. 
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Historically,  temperature  predictions  have  been  relatively  insensitive 
to  diffusion  formulations  when  compared  with  other  water  quality  param¬ 
eters  because  temperatures  are  physically  constrained  between  0°  and 
35°  C  by  heat  transfer  at  the  air/water  interface,  and  gradients  are 
usually  no  more  than  a  few  degrees  Celsius  per  metre.  In  contrast,  the 
magnitudes  and  gradients  of  other  water  quality  parameters  may  vary  over 
several  orders  of  magnitude  and  may  not  be  as  physically  constrained  to 
a  specific  range.  Model  simulations  of  these  parameters  can  therefore 
be  sensitive  to  changes  in  diffusion  coefficients  (Thornton  et  al. 

1979). 

Third,  the  algorithm  must  be  able  to  predict  changes  in  the  mixing 
regime  resulting  from  changes  in  hydrometeorological  conditions  and 
project  operation.  As  previously  stated  (Section  3.2),  mixing  in  res¬ 
ervoirs  results  from  hydrometeorological  forcing.  Superimposed  on  the 
seasonal  variation  in  hydrometeorological  forcing  are  short-term  fluct¬ 
uations.  The  mixing  regime  responds  to  both  of  these  forcing  scales. 
Different  mixing  regimes  result  depending  on  the  timing  of  the  forcing 
mechanism  and  the  formation  of  stratification.  A  change  in  project 
operation  also  results  in  a  change  in  the  mixing  regime.  For  example, 
increasing  the  pool  level  may  change  the  mixing  regime  from  an 
advective-dominated  to  a  buoyancy-dominated  regime.  Lowering  the  outlet 
depth  from  the  surface  to  the  bottom  increases  hypolimnetic  mixing.  The 
algorithm  must  be  able  to  simulate  these  changes  in  mixing  regime,  with¬ 
out  coefficient  changes,  in  order  to  be  used  for  evaluation  of  manage¬ 
ment  alternatives. 


5.2  Evaluation  Procedures 


Two  approaches  were  used  in  this  study  to  evaluate  alternative  mix¬ 
ing  algorithms: 

a.  Investigation  of  physical  basis. 

b.  Applications  to  different  reservoirs. 

Although  our  knowledge  of  reservoir  mixing  is  more  qualitative  than 
quantitative,  the  potential  sources  of  TKE  for  mixing  and  the  basic 
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functional  relationships  are  well  known  (Section  3.2).  Based  on  the 
review  of  reservoir  mixing  processes  in  Section  3.2,  the  important 
sources  of  TKE  were  identified  to  be  the  wind,  inflows,  and  outflows 
(including  location)  and  solar  and  atmospheric  heating  and  cooling 
(i.e.,  convective  mixing).  It  is  also  shown  that  mixing  increases  non- 
linearly  with  increases  in  wind  speed  and  inflow  rate.  Improvements  to 
existing  formulations  and/or  alternative  formulations  were  therefore 
first  evaluated  on  how  realistically  they  approximated  the  physical 
phenomena  to  be  modeled.  For  example,  several  of  the  formulations  for 
eddy  diffusion  coefficients  (i.e..  Equation  32)  were  discarded  since 
they  did  not  accurately  portray  mixing  processes  and/or  did  not  support 
field  observations  (Section  4.3.2). 

After  the  proposed  improvements  or  alternative  formulations  were 
determined  to  be  physically  viable,  they  were  evaluated  by  applications 
to  a  number  of  different  reservoirs  and  hydrometeorological  conditions. 
When  comparing  model  simulations  with  field  data,  it  is  important  to 
consider  the  previous  discussion  on  data  interpretation  (Section  2.2.5), 
especially  the  uncertainty  associated  with  field  measurements.  For 
example,  in  Figure  47  there  is  up  to  a  3°  C  difference  in  metalimnetic 
temperature  during  a  3-day  period.  Since  this  was  a  period  of  calm 
weather  and  low,  steady  inflows  and  outflows,  there  Is  little  reason  to 
expect  seiching  and  Internal  waves  that  could  cause  larger  temperature 
deviations  in  the  metalimnion.  The  difference  in  temperatures  is  prob¬ 
ably  representative  of  field  measurement  errors  and  should  be  considered 
when  comparing  field  data  with  model  predictions. 

Based  on  the  review  of  the  influence  of  mixing  on  reservoir  water 
quality  (Section  3.3),  it  was  determined  that  the  mixing  algorithm 
should  be  able  to  accurately  predict  the  onset  of  stratification,  mixed- 
layer  depth  dynamics,  metalimnetic  gradient,  variable  mixing  In  the 
hypolimnion,  fall  overturn,  inverse  stratification  during  winter  months, 
and  effects  of  project  operation  to  adequately  address  reservoir  water 
quality  problems. 

An  accurate  representation  of  the  onset  of  stratification  includes 
not  only  the  timing  but  also  the  depth  and  gradient.  As  explained  In 
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Figure  47.  Daily  variation  in  temperature  profiles, 

DeGray  Lake,  Arkansas 

Section  2.2  and  shown  in  Figure  48,  reservoir  stratification  starts  at 
the  bottom  of  a  lake  and  moves  upward,  not  from  the  water  surface  down¬ 
ward  as  some  older  diffusion  models  predict.  In  many  shallow 
reservoirs,  it  is  possible  for  stratification  to  form  and  break  up  sev¬ 
eral  times  prior  to  the  permanent  summer  stratification  forming. 

Daily  variations  in  mixed-layer  depths  result  from  synoptic  varia¬ 
tions  in  hydrometeorological  forcing  as  well  as  diel  variations  (Sec¬ 
tion  2.2.3).  Since  the  one-dimensional  assumption  limits  the 
computational  time  interval  to  1  day  (Section  4.2),  it  is  not  possible 
to  simulate  diel  variation  in  mixed-layer  dynamics.  Daily  variations  in 
mixed-layer  depths  and  temperatures  resulting  from  synoptic  variations 
in  hydrometeorological  forcing  are,  however,  significant.  As  indicated 
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in  Section  2.2.5,  it  is  sometimes  difficult  to  separate  diel  variations 
from  synoptic  variations  because  most  temperature  data  are  collected  at 
midday,  not  in  the  early  morning  hours. 

Since  density  gradients  inhibit  mixing  in  the  metalimnion,  the 
dynamics  of  the  metalimnion  are  sometimes  ignored.  An  accurate  repre¬ 
sentation  of  the  metalimnetic  gradient  requires  first  an  accurate  por¬ 
trayal  of  the  onset  of  the  stratification  and,  second,  the  inclusion  of 
mixing  from  synoptic  hydrometeorological  forcing.  It  is  important  to 
simulate  the  sharpening  of  the  metalimnetic  gradient  due  to  entrainment 
(Figure  49)  as  well  as  the  weakening  of  the  gradient  due  to  advection 
and  shear  (Figure  50) . 

Variable  mixing  in  the  hypollmnion  is  dependent  on  hydrometeorolog¬ 
ical  forcing  and  project  operation.  Although  little  is  known  about 
hypolimnetic  mixing,  it  is  known  that  mixing  levels  Increase  with  hydro¬ 
meteorological  forcing  and  use  of  bottom  gates.  Because  hypolimnetic 
temperatures  can  remain  relatively  constant,  it  is  difficult  to  deter¬ 
mine  the  significance  of  hypolimnetic  mixing  using  only  temperature 
data. 

With  fall  overturn,  complete  vertical  mixing  returns  to  the  lake. 
Consideration  must  be  given  to  both  the  water  temperature  (i.e.,  heat 
transfer,  Section  3.2.1)  and  the  magnitude  of  entrainment  (Figure  51). 
Accurate  simulation  of  this  process  may  require  the  inclusion  of  pene¬ 
trative  convection  (Section  2.4.6),  Because  of  the  triangular  longi¬ 
tudinal  profile  of  most  reservoirs  (Figure  52) ,  fall  overturn  appears  to 
start  at  the  upstream  end  of  the  reservoir  and  move  downstream  toward 
the  dam. 

Many  CE  reservoirs  are  located  in  cold  climates  and  experience  ice 
and  snow  covers  that  isolate  the  lakes  from  exchanges  across  the  air / 
water  interface.  The  lakes  then  become  inversely  stratified,  and  mixing 
is  dominated  by  inflow  and  outflow  processes.  The  actual  formation  of 
complete  ice  cover  may  take  as  little  time  as  a  few  hours  during  a  cold, 
calm  night  in  small  reservoirs,  to  several  weeks  in  large  reservoirs. 
Some  reservoirs  may  never  experience  complete  ice  cover.  To  accurately 
simulate  mixing  during  the  period  of  Ice  formation,  a  measure  of  the 
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Figure  49.  Sharpening  of  the  metalimnetic 
gradient  duo  to  entrainment,  DeGray  Lake, 
Arkansas 
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Figure  50.  Weakening  or  diffusion  of  the  metalimnetic 
gradient  due  to  shear  (inflow),  DeGray  Lake,  Arkansas 


Figure  51.  Fall  overturn  in  DeGray  Lake,  Arkansas 
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Figure  52.  Triangular  shape  of  a  reservoir 
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size  of  ice-free  water  is  required,  necessitating  the  use  of  partial 
area  ice  cover  algorithm  (Schultz  International,  Ltd.  1984). 

As  previously  described  (Section  3.2.6),  project  operation  can 
alter  the  mixing  regime  in  a  reservoir.  In  addition  to  considering  the 
nonlinear  effects  of  localized  TKE  input  and  its  interactions  with  den¬ 
sity  stratification,  an  accurate  representation  of  the  mixing  regime 
must  consider  other  algorithms  that  compute  effects  of  project  operation 
such  as  selective  withdrawal  zones. 


5.3  Historical  Development 

5.3.1  WQRRS 

Initial  work  on  the  mixing  algorithm  began  with  the  1974  version  of 
the  Water  Quality  for  River-Reservoir  Systems  model  (WQRRS)  (Hydrologic 
Engineering  Center  1974.)  WQRRS  was  the  basis  from  which  CE-QUAL-R1 
evolved  (Environmental  Laboratory  1982).  WQRRS  used  the  eddy  diffusion 
coefficient  formulation  shown  in  Figure  53.  Eddy  diffusion  coefficients 
were  constant  and  equal  in  both  the  epilimnion  and  hypolimnirn.  Meta- 
limnetic  coefficients  were  decreased  based  on  the  density  gradient  and  a 
calibration  coefficient,  A3.  The  metalimnion  was  distinguished  form  the 
epilimnion  and  hypolimnion  using  a  user-specified  stability  criterion, 
GSWH.  The  major  limitations  of  this  formulation  were  the  assump¬ 
tions  of  equal,  constant  eddy  diffusion  coefficients  in  the  epilimnion 
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Figure  53.  WQRRS  diffusion  coefficient 
formulations 
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and  hypolimnion  and  no  consideration  of  the  wind  as  a  mixing  mechanism 
(Figure  53b) . 

WQRRS  did  contain  an  alternative  formulation  that  considered  the 
wind  (Figure  53c) ,  but  in  this  formulation  the  dependence  on  the  wind 
speed  was  linear  and  therefore  not  physically  correct.  This  formulation 
also  did  not  consider  the  reduced  mixing  in  the  metalimnion  resulting 
from  density  gradients.  As  discussed  in  Section  4,3,2,  the  need  for 
reduced  mixing  in  the  metalimnion  is  critical,  and  models  that  do  not 
reduce  metalimnetic  eddy  diffusion  coefficients  with  increasing  density 
gradients  characteristically  predict  a  metalimnetic  gradient  that  is  too 
weak  even  if  the  simulations  are  started  with  a  measured  temperature 
profile  after  stratification  forms  (e.g..  Figure  54). 


Figure  54.  Diffused  metalimnetic  gradient  resulting 
from  using  a  constant  diffusion  coefficient 

Both  of  these  formulations  received  numerous  applications  and,  for 
the  most  part,  did  an  adequate  job  of  representing  the  thermal  regime  in 
reservoirs.  In  many  simulations  (e.g,,  Figure  55),  hypolimnetic  tem¬ 
peratures  were  too  low  in  the  spring  and  too  high  in  late  summer,  indi¬ 
cating  Insufficient  mixing  during  spring  turnover  and  too  much 
hypolimnetic  mixing  during  the  summer  stratified  period. 

In  order  to  eliminate  these  deficiencies  and  Improve  model  simu¬ 
lations,  the  formulation  in  WQRRS  was  modified  to  include  the  effects  of 
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Figure  55,  General  characteristics  of  predictions 
using  WQRRS,  Stockton  Lake,  Missouri  (to  convert 
feet  to  metres,  multiply  by  0.3048;  to  convert 
degrees  Fahrenheit  to  degrees  Celsius,  use  the 
following  formulas  C  ■  (5/9) (F  -  32)) 

the  wind  and  to  separate  mixing  in  the  epilimnion  from  mixing  in  the 
hypolimnion  (Figure  56),  The  formulation  retained  the  stability 
criteria  and  metalimnetic  formulation  from  WQRRS,  The  hypolimnetic 
mixing  coefficient  was  a  separate  user-specified  diffusion  coefficient 
that  was  considered  constant  throughout  the  simulation  period;  the  epi- 
limnetic  mixing  coefficient  was  dependent  on  the  wind  shear  stress 
(i.e.,  wind  speed  squared)  and  either  decreased  linearly  with  depth  to 
the  hypolimnetic  diffusion  coefficient  value  (Figure  56a)  or  was  con¬ 
stant  with  depth  (Figure  56b).  Simulations  of  a  number  of  different 
reservoirs  indicated  the  linear  decreasing  formulation  was  preferable. 
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PARAMETER  DIFFUSIVITY 


a.  Linear  with  depth 


b.  Constant  with  depth 

Figure  56,  Modified  WQRRS  diffusion  coefficient  formulation 
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As  with  WQRRS,  the  stability  criterion  GSWH  was  used  to  define  the 
density-dependent  metalimnion. 

This  eddy  diffusion  coefficient  formulation  considered  similar  phys¬ 
ical  processes  to  those  in  the  Richardson  number  formulations  discussed 
in  Section  4.3.2  but  had  more  coefficients  to  calibrate  to  match  field 
data.  Stafford  (1978)  showed  that  although  he  could  match  the  tempera¬ 
ture  structures  of  Lakes  Greeson,  Arkansas,  and  Sardis,  Mississippi,  for 
any  single  year,  it  was  not  possible  to  match  temperature  structures  for 
5  consecutive  years  using  the  same  calibration  coefficients.  Selected 
calibration  profiles  from  Lake  Greeson  (from  1972)  are  compared  with 
verification  simulations  from  1974  in  Figure  57.  It.  is  readily  apparent 
from  Figure  57  that  the  formulation  was  either  not  correctly  calibrated 
or  not  capable  of  simulating  different  hydrometeorological  conditions 
without  recalibration.  Simulations  using  other  Richardson  number-based 
formulations  produced  similar  results,  indicating  that  a  Richardson 
number-type  formulation  for  the  diffusion  coefficients,  alone,  was  not 
adequate  to  accurately  simulate  changes  in  a  lake's  thermal  structure 
resulting  from  hydrometeorological  conditions. 

5.3.2  Variable  Layer  Formulation 

One  of  the  problems  associated  with  the  1974  Lake  Greeson  predic¬ 
tions  (Figure  57)  was  excessive  metalimnetic  mixing.  This  mixing 
resulted,  in  part,  from  the  vertical  advective  term  (term  (2))  in  Equa¬ 
tion  29  for  two  reasons.  First,  reduction  of  the  mixing  coefficients  to 
produce  molecular  diffusion  in  this  region  did  not  improve  model  predic¬ 
tions.  Second,  the  outlet  was  located  in  the  metalimnion  at  the  same 
elevation  where  excessive  mixing  was  observed. 

To  reduce  this  numerical  mixing,  the  variable  layer  concept  or 
Lagrangian  approach  of  Imberger  et  al.  (1978)  was  incorporated  into  the 
model.  With  this  approach  there  is  no  vertical  advection  term  and  the 
layers  are  allowed  to  expand  or  contract  with  layer  inflows  and  out¬ 
flows.  A  simple  example  illustrating  the  differences  between  the  fixed 
layer  model  and  variable  layer  model  is  shown  in  Figure  58.  With  the 
fixed  layer  model,  the  net  result  of  the  mass  Inflow  to  the  top  layer 
and  withdrawal  from  the  bottom  layer  is  an  increase  in  bottom  layer 
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Figure  58.  Schematic  comparing  fixed  and  variable 
layer  formulations 

concentration  (0.5C  to  0.6C)  and  decrease  in  the  vertical  concentration 
gradient  between  the  two  layers.  The  volumes  or  sizes  of  the  two  layers 
remain  constant.  In  contrast,  with  the  variable  layer  model,  the  top 
layer  increases  in  size  and  the  bottom  layer  decreases  in  size  while  the 
layer  concentrations  and  the  gradient  remain  unchanged. 

A  comparison  of  model  predictions  using  the  variable  and  fixed 
layer  formulations  is  shown  in  Figure  59.  As  expected,  there  was  less 
metalimnetic  mixing  with  the  variable  layer  formulation. 

5.3.3  Mixed-Layer  Dynamics 

In  addition  to  problems  associated  with  simulating  consecutive 
years  on  Lake  Greeson,  the  modified  WQRRS  formulation  also  did  not  accu¬ 
rately  portray  the  onset  of  stratification  and  mixed-layer  dynamics. 
Since  Ford  (1976)  showed  that  a  simple  TKE  formulation  (Section  4.4)  was 
able  to  accurately  reproduce  the  onset  of  stratification  and  mixed-layer 
dynamics  in  small  natural  lakes  and  Hurley-Octavlo ,  Jirka,  and  Harleman 
(1977)  and  Imberger  et  al.  (1978)  were  able  to  simulate  the  thermal 
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structure  of  larger  reservoirs  using  similar  formulations,  a  TKE 
formulation  or  mixed-layer  formulation  was  incorporated  into  the  mixing 
algorithm. 

Initially,  the  mixed  layer  formulation  followed  the  approach  of 
Ford  (1976)  and  Hurley-Octavio,  Jirka,  and  Harleman  (1977)  and  did  not 
consider  dissipation  of  TKE  with  depth.  Although  the  simulations  of  the 
onset  of  stratification  improved  (Ford  et  al.  1981),  there  was  too  much 
entrainment  during  periods  of  weak,  stratification  (i.e.,  the  predicted 
mixed- layer  depths  were  too  deep)  (Figure  60).  Imberger  et  al.  (1978) 
avoided  this  problem  by  dissipating  the  TKE  with  depth  using  a  cubic 
function.  There  was,  however,  no  physical  basis  for  the  function.  Not¬ 
ing  the  same  problem  with  the  Hurley-Octavio  model,  Bloss  and  Harleman 
(1980)  developed  an  efficiency  function  for  entrainment  based  on  the 
Richardson  number  (Section  4.4).  This  function  significantly  improved 


Figure  60.  Comparison  of  simulation  results  from  a  diffusion  and 
TKE  model,  Lake  Anna,  Virginia  (after  Hurley-Octavio,  Jirka,  and 

Harleman  1977) 
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model  simulations  (Figure  61).  Johnson  and  Ford  (1981)  were  able  to 
simulate  6  years  of  data  at  DeGray  Lake,  Arkansas,  and  5  years  at  Lake 
Greeson,  Arkansas  (including  the  years  Stafford  (1978)  had  problems 
with).  Figures  62  and  63  compare  Stafford's  results  using  a  diffusion 
model  with  Johnson  and  Ford's  results  using  the  TKE  formulation  in  a 
diffusion  model. 

In  the  1975  simulations  (Figure  62),  the  onset  of  stratification 
(Julian  Day  (JD)  91),  mixed-layer  depths,  metalimnetic  gradients,  and 
hypolimnetic  temperatures  (JD  167,  251,  and  321)  were  all  more  accu¬ 
rately  simulated  with  the  mixed-layer  formulation.  The  1973  simulations 
(Figure  63)  were  also  significantly  better  with  respect  to  these  param¬ 
eters.  As  shown  in  Figure  64,  the  mixed-layer  depth  or  TKE  formulation 
was  able  to  accurately  predict  differences  in  metalimnetic  gradients 
resulting  from  hydrometeorological  forcing.  As  previously  indicated 


Temperature  ,*C 


Figure  61.  Comparison  of  simulation  results  with  and  without 
dissipation  function,  Lake  Anna,  Virginia  (after  Bloss  and 

Harlemann  1980) 
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Figure  64.  Simulation  o:  metallmnetic  gradient,  Lake  Greeson,  Arkansas 

(Section  5.3.1),  pure  difrusLon  models  based  on  Richardson  number  for¬ 
mulations  were  not  able  to  reproduce  these  differences. 

5.3.4  Diffusion  Coefficient  Formulation 

The  incorporation  of  the  variable  layer  concept  and  mixed-layer 
dynamics  into  the  model  reduced  numerical  dispersion  (Section  5.3.2), 
improved  mixed-layer  depth  predictions  (Section  5.3.3),  and  identified 
the  need  to  modify  the  diffusion  coefficient  formulacion.  Model  simu¬ 
lations  using  the  modified  formulation  indicated  that  mixing  from 
inflows  and  outflows  was  not  being  adequately  considered  and  that  the 
hypolimnetic  mixing  coefficient  should  be  variable  and  dependent  on  the 
wind  speed. 

Analysis  of  storm  event  data  from  DeGray  Lake,  Arkansas,  clearly 
showed  that  mixing  from  inflows  diffuses  or  weakens  the  metallmnetic. 
temperature  gradient  (Figure  50).  In  the  fixed  layer  formulation, 
mixing  from  inflows  and  outflows  was  introduced  through  the  numerical 
dispersion  introduced  from  the  vertical  advection.  Elimination  of  this 
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dispersion  with  the  incorporation  of  the  variable  layer  concept  also 
eliminated  much  of  the  mixing  associated  with  inflows  and  outflows  in 
the  model.  Since  Ozmidov  (1965)  related  the  eddy  diffusion  coefficient 
to  the  dissipation  rate  of  TKE  (Equation  44) ,  mixing  from  inflows  and 
outflows  was  incorporated  into  the  eddy  diffusion  coefficient  by  com¬ 
puting  the  TKE  input  associated  with  advection  in  a  horizontal  layer 
i  : 


TKE  =  -Jr  p  Q,  V2.  (50) 
a  2  w  i  i 

where 

3 

=  flow  rate  in  layer  i,  m  /sec 
Ui  =  flow  velocity  in  layer  i,  m/sec 

Equation  44  was  used  to  incorporate  variable  mixing  into  the  hypo- 
limnion  of  a  lake.  Since  wind-generated  TKE  can  indirectly  enter  the 
metalimnior  and  hypolimnion  via  seiche  motion  and  breaking  of  internal 
waves,  the  wind-generated  TKE  was  assumed  to  be  dissipated  throughout 
the  entire  reservoir. 

These  formulations  were  incorporated  in  CE-QUAL-Rl  and  used  by 
Johnson  and  Ford  (1983)  to  simulate  the  thermal  structures  of  Lakes 
DeGrav  and  Greeson, 


5.4  Recommended  Algorithm 


Based  on  an  investigation  of  the  physical  basis  for  a  number  of 
algorithms  and  applications  to  reservoirs  representing  different  regions 
of  the  country,  different  hydrometeorological  conditions,  and  different 
operating  conditions,  it  is  recommended  that  the  generalized,  one- 
dimensional  model  capable  of  simulating  changes  in  the  mixing  regime 
include : 

a.  Variable  layer  formulation. 
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b.  Mixed-layer  dynamics  using  the  efficiency  function  of  Bloss  and 
Harleman  (1980)  and  considering  penetrative  convection. 

c.  Variable  diffusion  coefficient  that  depends  on  wind,  inflows, 
outflows,  and  density  differences. 

5.4.1  Variable  Layer  Formulation 

A  Lagrangian  or  variable  layer  formulation  is  recommended  to  reduce 
numerical  diffusion  and  to  allow  calibration  of  the  mixing  resulting 
from  inflows  and  outflows.  Although  the  details  and  advantages  of  the 
variable  formulation  using  power  functions  for  the  area  A(z) 


$ 


A(z)  »  a^  Z  exp  ( a 


and  the  volume  V(z) 


al 

V(z)  *  - — 7.  exp  (a2  +  1) 


are  fully  described  in  the  CE-QUAL-Rl  User's  Manual  (Environmental 
Laboratory  1982),  experience  has  shown  that  these  simple  power  functions 
do  not  always  provide  a  sufficiently  accurate  representation  of  reser¬ 
voir  geometry  for  water  quality  simulations.  It  is  therefore  recom¬ 
mended  that  a  polynomial  function  or  a  lookup  table  of  actual  areas  and 
volumes  be  used  in  future  versions  of  CE-QUAL-Rl. 

5.4.2  Mixed-Layer  Dynamics 

A  mixed  layer  of  TKE  formulation  is  recommended  to  accurately  simu¬ 
late  the  onset  of  stratification  and  mixed-layer  dynamics.  This  type  of 
formulation  ensures  that  stratification  starts  at  the  bottom  of  the  lake 
and  moves  up,  and  that  the  dynamics  of  hydrometeorological  forcing  deter¬ 
mine  the  shape  of  the  metalimnetic  temperature  gt  dient.  An  accurate 
representation  of  mixed-layer  dynamics  during  periods  of  cooling 
requires  consideration  of  penetrative  convection. 

Wind .  The  TKE  available  for  possible  entrainment  from  the  wind 
shear  can  be  estimated  by 
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(53) 


TKE 

w 


wJ 


T 

S 


At  dA 


where 

2  2 

TKE  »  wind  shear  turbulent  kinetic  energy,  (kg-m  )/sec 
w  2 

A  “  water  surface  area,  m 

p 

c -  empirical  coefficient 

w*  “  shear  velocity  of  water,  m/sec 

2 

x  »  shear  stress  at  the  air/water  interface,  kg/(m-sec  ) 
s 

At  ■  time  step,  sec 


The  shear  velocity  w^  is  defined  by 


w 


I 


in  which  p 

w 


water  density,  kg/m 


3 


The  shear  stress  at  the  air /water  interface  is  given  by: 


T 

S 


n 

'  a 


where 

p  -  air  density  (1.177  kg/m^) 

a 

®  dimensionless  drag  coefficient 
W  =  wind  speed,  m/sec 


(54) 


(55) 


The  drag  coefficient  is  taken  from  Safaie  (1978): 

C,  -  0.00052  W0’44  (56) 

d 

In  many  lakes,  the  wind  stress  does  not  act  on  the  entire  surface  area 
because  of  sheltering  effects  from  the  surrounding  terrain.  The  TKE^ 
(Equation  53)  Is  therefore  modified  by  a  site-dependent  sheltering 
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coefficient  that  is  the  ratio  of  the  watsr  surface  area  exposed  to  the 
wind  to  the  total  water  surface  area*  The  coefficient  has  a  maximum 
value  of  L  when  sheltering  is  insignificant.  It  may  also  be  necessary 
in  some  regions  to  modify  Equation  56  for  the  drag  coefficient  to  con¬ 
sider  nonneutral  atmospheric  conditions. 

Penetrative  convection.  TKE  produced  by  convection  currents  during 
periods  of  cooling  is  assumed  to  be  proportional  to  the  net  heat  flux 
,  when  is  negative.  The  energy  available  for  entrainment  from 

overturn  convection  can  be  estimated  by: 

TKE  -  -c  H  A  h  g  a  At/c  (57) 

c  c  n  s  p 


where 

2  2 

TKE^  **  turbulent  kinetic  energy,  (kg-m  )/sec 

c  *  empirical  calibration  coefficient 

C  2 
*  net  heat  flux  across  the  air/water  interface,  W/m 

h  *  depth  of  mixed  layer,  m 

2 

g  •»  acceleration  due  to  gravity „  m/sec 
a  «*  coefficient  of  thermal  expansion  for  water,  per  °C 
0^  specific  heat  of  water,  J/(kg-°C) 

In  Equation  57,  TKE^  equals  0  when  is  positive.  The  total  TKE 

available  for  entrainment  is: 

TKE  -  TKE  +  TKE  (58) 

w  c 

Mixing  ef ficlenclee.  Because  mixing  processes  are  dissipative  and 
not  efficient  and  because  different  processes  dominate  at  different 
tim^s,  Equation  58  must  be  modified.  Based  on  parameterization  of  the 
TKE  balance  at  a  density  interface,  Bloss  and  Harleman  (1980)  determined 
a  Richardson  number  function  f(R^)  to  modify  the  TKE: 

29.46  -  YRT 

«V  ■  °-057  Ri  -too  (59) 
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where 


2 

R.  ■  Richardson  number  re  hgAp/p  w. 
i  w  * 

Ap  ■  density  difference  across  the  interface 

The  total  TKE  available  for  entrainment  is 

TKE  =  TKE  *  F(Ri)  (60) 

Entrainment.  Prior  to  calculating  entrainment,  a  temporary  tem¬ 
perature  (density)  structure  is  computed  for  the  computation  interval. 
This  structure  considers  internal  absorption  of  solar  radiation,  net 
heat  transfer  at  the  air /water  interface,  inflows,  and  outflows.  It 
does  not  consider  mixing  between  layers. 

The  work  W  required  to  entrain  or  lift  the  mass  ApAV  from  its 
position  immediately  below  the  mixed  layer  to  the  center  of  mass  of  the 
mixed  layer  is 

WT  -  Ap  AV  g  (h  -  h  )  (61) 

where 

2  2 

«  entrainment  work,  (kg-m  ) /sec 
Ap  -  density  difference  between  the  mixed  layer  and 

3 

underlying  layer,  kg/m 

3 

AV  <=  incremental  volume  to  be  entrained,  m 

2 

g  *  acceleration  due  to  gravity,  m/sec 
h  *  depth  of  the  center  of  mass  of  the  mixed  layer,  m 

The  depth  of  the  mixed  layer  is  calculated  after  comparing  Equations  60 
and  61.  If  the  TKE  is  larger  than  W  ,  entrainment  occurs  and  h 

if 

increases.  Entrainment  continues  until  TKE  is  no  longer  larger  than 
WL  ' 

As  indicated  in  Section  4.4.2,  this  type  of  formulation  always 
sharpens  the  metalimnetic  gradient.  Diffusion  or  weakening  of  the  gra¬ 
dient  is  discussed  in  the  next  section. 
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5.4,3  Diffusion  Coefficient  Formulation 

Following  the  work  of  Ozmidov  (1965),  the  recommended  form  of  the 
diffusion  coefficient  for  stratified  conditions  is 


K 

z 


1  N2 


(62) 


where 


K 


N 


'1 

e 

,2 


2 

global  vertical  diffusion  coefficient,  m  /sec 
dimensionless  calibration  coefficient 


2  3 

local  dissipation  rate  of  TKE,  m  /sec 


g  9p 
p  3z 


buoyancy  frequency,  sec 


-2 


This  equation  assumes  a  local  steady  state  or  equilibrium  such  that  the 
rate  of  input  of  TKE  at  larger  scales  is  in  balance  with  the  rate  of 
dissipation  at  smaller  scales.  Equation  62  has  been  theoretically  jus¬ 
tified  by  Weinstock  (1978),  who  assumed  mixing  was  done  in  the  inertial 
subrange.  These  assumptions  are  acceptable  since  turbulence  is  highly 
dissipative.  Values  for  vary  from  0,8  (Weinstock  1978)  to  0.25 

(Linden  1979,  McEwan  1980,  Oakey  1982).  Equation  62  states  that  mixing 
or  the  diffusion  coefficient  increases  with  increasing  energy  input  and 
decreases  with  increasing  density  gradient. 

The  computation  of  the  local  dissipation  rate  needs  to  consider  the 
TKE  Inputs  from  the  wind,  inflows,  and  outflows  separately,  since  the 
inputs  from  the  inflows  and  outflows  are  not  spread  uniformly  over  the 
entire  lake  as  the  wind  is  but  are  restricted  to  the  zone  of  inflow  or 
outflow  (Imberger  and  Patterson  1981).  This  finding  was  also  verified 
in  simulations  of  DeGray  Lake, 

The  total  TKE  input  rate  from  the  wind  is  given  by  Equation  53. 

The  rate  of  energy  dissipation  for  wind  per  unit  mass  is 
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a 


(63) 


P„  v*  A 


w 


P  V 
w 


D 

m 


where 

V  *  volume, 

D  ■  mean  depth,  m 
m 

This  relationship  is  similar  to  the  rate  of  dissipation  in  an  unstrati¬ 
fied  uniform  stress  layer  where  shear-generated  turbulence  cascades  to 
dissipation  scales  (e.g.,  Turner  1973).  This  is 


e 


3U  - 

-  Be 

dz 


3 

w* 

kz 


(64) 


where 

u  u  -  Reynolds  stress,  m^/sec^ 

X  z 

3U 

—  *  velocity  gradient,  (m/sec) /m 

d  Z 

k  ■«  von  Kantian' s  constant  (~0.4) 

This  relationship  has  been  found  to  be  valid  in  the  surface  layer  of  a 
reservoir  (Dillon  et  al.  1981)  and  in  the  mixed  layer  off  the  Nova 
Scotia  continental  shelf  (Oakey  and  Elliott  1980,  1982). 

The  TKE  input  from  the  inflow  and  outflow  is  given  by  Equation  50. 
The  dissipation  rate  is  therefore 


e 

a 


TKEa  1  K  Ui\ 
\  Vi  “  5  \  Vi  / 


where  *=  volume  of  layer  i,  m  . 


(65) 


With  the  substitution  of  Equations  63  and  65  into  Equation  62,  the  dif¬ 
fusion  coefficient  becomes 
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(66) 


K 


z 


c 


w 


£ 

w 


+  c 

a 


N 


2 


e 

a 


where 

c  =  calibration  coefficient  for  wind  mixing,  dimensionless 
w 

*  calibration  coefficient  for  advective  mixing,  dimensionless 

There  are  two  problems  with  Equation  66.  First,  the  diffusion  coeffi- 
cient  becomes  infinitely  large  when  the  density  gradient  (i.e.,  N  )  goes 
to  zero  (i.e.,  a  well-mixed  condition).  Second,  the  diffusion 
coefficient  is  dependent  on  the  density  gradient  to  the  -1  power,  but 
field  data  indicate  this  power  varies  between  -1/2  and  -2  (Sec¬ 
tion  4.3.2).  These  deficiencies  were  corrected  by  defining  the  dimen¬ 
sionless  stability  factor 


<.  ,  1L  la 

f  Ap  3z 


(67) 


where 

H  ■  depth  of  lake,  m 

3 

Ap  *  density  differential  between  bottom  and  surface  waters,  kg/m 

For  well-mixed  conditions,  ■=  1  .  Incorporation  of  into  Equa¬ 

tion  66  requires  the  addition  of  a  time  constant  to  keep  Equation  66 
dimensionally  correct.  This  time  constant  was  assumed  to  be  the  compu¬ 
tation  interval.  The  formulation  for  the  diffusion  coefficient  is 
therefore 


K 

z 


(sf)n 


(68) 


where  n  =  calibration  coefficient  . 
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This  formulation  for  the  diffusion  coefficient  has  several  desir¬ 
able  features: 

a.  Mixing  increases  with  energy  input. 

b.  Mixing  resulting  from  the  wind  and  flow  can  be  calibrated 
separately. 

c.  The  formulation  is  density  dependent. 

d.  The  density  dependency  can  be  varied. 

5.5  Applications  and  Verification 


During  the  period  1976  to  1983,  the  recommended  algorithm  and  its 

predecessors  were  used  to  simulate  the  thermal  structures  of  over 

15  reservoirs  and  lakes  and  numerous  proposed  reservoirs  (i.e.,  over 

40  data  sets).  A  partial  listing  is  given  in  Table  3.  The  reservoirs 

and  lakes  are  from  seven  states  and  one  province  in  Canada  and  represent 

all  types  of  hydrometeorological  conditions.  The  lakes  vary  in  size 

from  F.  E.  Walter  Reservoir,  Pennsylvania,  which  has  a  volume  of  2.47  x 
6  3 

10  m  and  mean  annual  hydraulic  residence  time  of  1.7  days  (Ford, 

Thornton,  and  Norton  1983),  to  Williston  Reservoir  which  has  a  volume  of 
10  3 

5.7  x  10  m  and  a  mean  annual  hydraulic  residence  time  of  several 
years  (Schultz  International,  Ltd.  1984). 

The  applications  clearly  illustrate  that  the  recommended  one¬ 
dimensional  algorithm  is: 

a.  Sufficiently  general  to  simulate  the  thermal  structures  of  all 
types  and  sizes  of  CE  reservoirs. 

b.  Capable  of  simulating  differences  in  a  lake's  thermal  structure 
resulting  from  hydrometeorological  conditions. 

c.  Capable  of  simulating  differences  in  a  lake's  thermal  structure 
resulting  from  changes  in  operation  (i.e.,  different  withdrawal 
depth,  different  rule  curve). 

d.  Capable  of  simulating  the  thermal  structure  of  a  proposed  res¬ 
ervoir  after  being  calibrated  on  a  morphometrically  similar, 
existing  reservoir. 

To  illustrate  the  last  point,  the  model  was  calibrated  on  1979  data  from 
DeGray  Lake  using  bottom  withdrawal  and  then  applied  to  Lakes  Greeson, 
Ouachita,  Hamilton,  and  Catherine.  These  four  lakes  are  located  within 


Table  3 

Summary  of  Model  Applications 


Reservoir 

Simulation 

Year(s) 

Purpose 

Beltzville  Lake, 
Pennsylvania 

1972 

Verification  of  1981  F.  E. 
Walter  calibration 

1981 

Verification  of  1981  F.  E. 
Walter  calibration 

Lake  Catherine,  Arkansas 

1982 

Verification  using  1979 

DeGray  calibration 

Lake  Calhoun,  Minnesota 

1974 

Calibration 

1975 

Verification 

C.  J.  Brown  Reservoir, 

Ohio 

1974 

Simulation  of  filling 

1975 

Calibration 

Lake  Coralville,  Iowa 

1966-67 

Calibration 

1969 

Calibration 

DeGray  Lake,  Arkansas 

1974-78 

Verification,  surface 
withdrawal 

1979 

Calibration,  bottom 
withdrawal 

1980 

Verification,  bottom 
withdrawal 

1982 

Verification,  bottom 
withdrawal 

F.  E.  Walter  Reservoir,. 
Pennsylvania 

1977 

Verification  of  1981  on 
shallow  pool,  16.8  m  deep 

1979 

Verification  of  1981  on 
shallow  pool,  16.8  m  deep 

1981 

Calibration  on  high  pool, 

44.8  m  deep 

Lake  Greeson,  Arkansas 

1972-76 

Verification  of  1979  DeGray 
calibration 

1982 

Verification  of  1979  DeGray 
calibration 

Lake  Hamilton,  Arkansas 

1982 

Verification  of  1979  DeGray 
calibration 

(Continued) 
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Table  3  (Concluded) 


Reservoir 

Simulation 
Year (s) 

Purpose 

McCarrons  Lake,  Minnesota 

1974 

Verification 

calibration 

of 

1974  Calhoun 

1975 

Verification 

calibration 

of 

1974  Calhoun 

Lake  Ouachita,  Arkansas 

1982 

Verification 

calibration 

of 

1979  DeGray 

Sardis  Lake,  Mississippi 

1966 

Verification 

calibration 

of 

1968 

1967 

Verification 

calibration 

of 

1968 

1968 

Calibration 

1969 

Verification 

calibration 

of 

1968 

1970 

Verification 

calibration 

of 

1968 

Lake  Shelbyville, 

Illinois 

1973 

Calibration 

1975 

Verification 

calibration 

of 

1973-77 

1977 

Calibration 

Williston  Lake,  British 
Columbia 

1976-77 

Verification 

calibration 

of 

1981-82 

1981-82 

Calibration 

50  km  of  DeGray  Lake  in  southwestern  Arkansas.  Their  major  morphometric 
characteristics  are  compared  in  Table  4.  Johnson  and  Ford  (1983)  dis¬ 
cuss  the  applications  to  several  years  of  data  for  Lake  Greeson  while 
FTN  (1983)  presents  the  results  for  Lakes  Ouachita,  Hamilton,  and 
Catherine. 
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Table  4 

Comparison  of  Major  Morphometric  Characteristics 


SECTION  6.  SUMMARY  AND  CONCLUSIONS 


6,1  Literature  Review  Findings 

The  review  of  reservoir  mixing  processes  indicated  the  observed 
thermal  structure  in  a  reservoir  results  from  the  cumulative  effects  of 
a  number  of  complex,  nonlinear,  interdependent  mixing  processes. 

Although  the  sources  of  TKE  for  mixing  are  limited  to  solar  and  atmo¬ 
spheric  heating,  wind,  inflows,  and  outflows  (including  outlet  loca¬ 
tion)  ,  quantitative  knowledge  of  specific  mixing  mechanisms  such  as  wind 
mixing,  penetrative  convection,  turbulent  diffusion,  and  Kelvin- 
Helmholtz  instabilities  is  limited.  The  TKE  budgets  for  the  entire  lake 
do,  however,  clearly  identify  the  physical  causes  of  observed  motion  and 
thermal  stratification.  They  also  indicate  which  source  of  energy  is 
controlling  mixing  at  a  particular  time.  It  is  common  for  wind  to  domi¬ 
nate  at  one  time  and  inflows  to  dominate  at  another  time. 

The  review  of  the  influence  of  mixing  on  reservoir  water  quality 
indicated  that  mixing  significantly  impacts  and  sometimes  controls  the 
observed  water  quality  by  controlling  horizontal  and  vertical  constit¬ 
uent  distributions  and  thereby  influencing  the  physical,  chemical,  and 
biological  regimes.  The  review  also  indicated  that  in  order  for  a  one¬ 
dimensional  model  to  accurately  simulate  the  effects  of  project  opera¬ 
tion  on  reservoir  water  quality,  the  mixing  algorithm  should  be  capable 
of  accurately  simulating  the: 

a.  Onset  of  stratification. 

b.  Daily  variations  in  mixed-layer  depths  and  dynamics  of  short¬ 
term  mixing  events. 

c.  Metalimnetic  gradient. 

d_.  Variable  mixing  in  the  hypo  limn  ion. 

e.  Fall  overturn. 

£.  Inverse  stratification  during  winter  months. 

Effects  of  project  operation. 
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The  review  of  one-dimensional  predictive  techniques  indicated  that 
the  recommended  mixing  algorithm  should  include  formulations  for  both 
diffusion  and  mixed-layer  dynamics  (i.e,,  TKE  models).  Inclusion  of 
only  one  type  of  formulation  severely  limits  the  applicability  of  the 
model  for  CE  reservoirs. 

6.2  Recommended  Algorithm 

Since  the  recommended  mixing  algorithm  is  intended  to  be  used  in  a 
generalized  one-dimensional  water  quality  model  (CE-QUAL-R1)  to  predict 
changes  in  reservoir  water  quality  resulting  from  changes  in  hydromete¬ 
orological  conditions  and  project  operation,  several  requirements  were 
considered  during  algorithm  development: 

a.  The  algorithm  had  to  be  generalized  with  respect  to  CE 
reservoirs  and  therefore  not  be  constrained  by  extensive 
morphometric  and  hydrometeorological  data  requirements  nor 
limited  in  the  mixing  processes  considered. 

b.  The  algorithm  is  to  be  used  in  a  one-dimensional  water  quality 
model  that  does  not  compute  the  internal  current  structure. 

c.  The  algorithm  must  include  all  major  mixing  processes  in  order 
to  predict  changes  in  the  mixing  regime  resulting  from  changes 
in  hydrometeorological  conditions  and  project  operation. 


Based  on  the  review  of  reservoir  mixing  processes,  an  investiga¬ 
tion  of  the  physical  basis  for  a  number  of  algorithms,  and  applications 
to  reservoirs  representing  different  regions  of  the  country,  different 
hydrometeorological  conditions,  and  different  operating  conditions,  the 
recommended  one-dimensional  algorithm  includes: 

a.  Variable  layer  formulation  (Section  5.4.1)  to  reduce  numerical 
dispersion  and  to  allow  direct  calibration  of  mixing  resulting 
from  inflows  and  outflows. 

b.  Mixed-layer  dynamics  using  the  TKE  available  for  possible 
entrainment  from  wind  shear  (Equation  53) »  TKE  produced  by 
convective  currents  during  periods  of  cooling  (Equation  57), 
and  the  efficiency  function  of  Bloss  and  Harleman  (1980) 
(Equation  59). 
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c.  Variable  diffusion  coefficient  that  depends  on  the  energy 

inputs  from  the  wind,  inflows,  and  outflows  and  on  the  density 
gradient  (Equation  68). 

The  recommended  algorithm  was  used  to  simulate  the  thermal  struc¬ 
tures  of  over  15  reservoirs  and  lakes  of  varying  geographical  location, 
size,  hydrometeorological  regime,  and  operation  configurations. 
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APPENDIX  A:  STRATIFICATION  COMPUTATIONS 

A.l  Introduction 


Most  lakes  and  reservoirs  stratify,  albeit  weakly  and  intermit¬ 
tently,  at  one  time  or  another.  The  major  factors  that  limit  strat¬ 
ification  are  the  lake  depth,  flowthrough  rate,  and  the  wind.  These 
factors  have  been  used  by  several  investigators  to  develop  criteria  for 
stratification  potential  (i.e.,  the  likelihood  that  a  particular  body  of 
water  will  stratify).  In  general,  lakes  with  depths  greater  than  10  m 
and  mean  annual  residence  times  greater  than  20  days  stratify.  These 
numbers  differ  slightly  from  Harleman  (1982),  who  states  that  lakes  with 
depths  greater  than  5  m  stratify  except  if  the  residence  time  is  less 
than  36  days.  In  addition  to  these  general  rules,  several  computations 
can  be  made  to  evaluate  the  stratification  potential  of  a  reservoir. 
Note,  however,  that  these  calculations  are  merely  indicators  and,  as 
shown  in  the  examples,  all  methods  do  not  always  yield  the  same 
conclusion. 


A. 2  Densimetric  Froude  Number 


A. 2.1  Description 

Norton,  Roesner,  and  Orlob  (1968)  proposed  a  more  scientifically 
based  stratification  criterion  in  the  form  of  the  densimetric  Froude 
number 


(Al) 


where 

F,  =  densimetric  Froude  number,  dimensionless 

d  2 
g  =  acceleration  of  gravity,  m/sec 

e  =  dimensionless  density  gradient,  10  ^/m 

L  *  reservoir  length,  m 


Al 


_  3 

Q  -  average  reservoir  discharge,  m  /sec 

D  -  reservoir  mean  depth,  in 

m  3 

V  =  reservoir  volume,  m 

An  >>  1/ir  indicates  a  well-mixed  system;  F^  «  1/tt  indicates  a 
strongly  stratified  system;  while  F^  -  1/tt  indicates  a  weakly  or 
intermittently  stratified  system. 

This  criterion  compares  the  destratifying  force  of  the  flowthrough 
with  the  stratifying  potential  of  the  assumed  density  gradient.  It  can 
be  simplified  and  rewritten  in  terms  of  the  residence  time  t 
(t  *=  volume/flow  rate) 
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319  L 
t  D 


r  m 


(A2) 


For  the  critical  F  of  1/tt  and  a  residence  time  of  20  days, 

d  -4 

L/D  »  1,724  or  the  bottom  slope  is  on  the  order  of  5.8  x  10  .  Since 

m 

this  slope  is  characteristic  of  many  CE  reservoirs,  the  critical 
residence  time  of  20  days  is  consistent  with  the  Norton,  Roesner,  and 
Orlob  (1968)  criterion. 

A. 2. 2  Examples 

DeGray  Lake  is  a  CE  impoundment  located  on  the  Caddo  River  in 
south-central  Arkansas.  It  has  the  following  characteristics: 


L  *  32,000  m 

D  “  14 . 8  m 
m 

V  «  7.91  x  10®  m3 

_  3 

Q  =  18.2m/ sec 

The  densimetric  Froude  number  is,  therefore. 


v  /IT  iL  =  f  1  (32,000)  (18. 

Vg  DmV  V  (9 . 8) ( 1U-6)  (14.8) (7.91  x 


2) _ 

108) 


0.02 
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Since  F  *  0.02  «  1/ir  ,  DeGray  Lake  should  be  a  strongly  stratified 
d 

system*  which  it  is  (see  Figure  6).  Its  maximum  depth  of  57  m  and  mean 
annual  residence  time  of  1.38  years  also  confirm  this  conclusion. 

F.  E.  Walter  Reservoir  is  a  small  CE  project  located  on  the  Lehigh 
River  in  the  Pocono  Mountains  of  northeastern  Pennsylvania.  It  has  the 
following  characteristics: 


L  -  2,700  m 

D  *  6.8  m 
m 

V  -  2.47  x  106  m3 
3 

Q  ■  17.6  m  /sec 
The  densimetric  Froude  number  is,  therefore, 


(2,700)  (17.6) 


(9.8) (10-6)  (6.8) (2.47  *  106) 


0.9 


indicating  a  weakly  or  intermittently  stratified  system.  In  contrast, 
its  maximum  depth  of  16.8  m  would  indicate  a  stratified  system  and  its 
mean  annual  residence  time  of  1.7  days  would  indicate  a  well-mixed 
system.  Since  F.  E.  Walter  Reservoir  intermittently  stratifies,  the 
densimetric  Froude  number  computation  gives  the  proper  result  while  the 
conflicting  maximum  depth  and  residence  time  criteria  compensate  one 
another. 


A. 3  Wind  Mixing  Depth 


A. 3 . 1  Description 

The  depth  of  the  lake  is  important  when  evaluating  stratification 
potential  since  the  effects  of  wind  mixing  and  the  penetration  of  solar 
radiation  are  depth  limited.  If  solar  radiation  does  not  penetrate  to 
the  bottom  of  the  lake,  wind  mixing  is  required  to  prevent  stratification 


A3 


from  forming.  Since  the  average  Secchi  disc  depth  of  CE  reservoirs  is 
1.1  m  (Thornton,  Nix,  and  Bragg  1980),  wind  mixing  is  required  in  most 
CE  reservoirs  to  mix  the  heat  to  depths  greater  than  2.2  m.  The 
importance  of  wind  mixing  can  be  evaluated  using  the  length  or  depth 
scale  (Sundaram  1973): 
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where 

w^  «  shear  velocity  of  wind  in  water,  m/sec 
8  =  empirical  coefficient,  dimensionless 

a  ■  volumetric  coefficient  of  thermal  expansion,  1.8  *  10  /° C 

2 

H  »  net  surface  heat  flux,  W/m 

n  3 

p  «  density  of  water,  kg/m 

0^  ■  specific  heat  of  water,  4,186  J/(kg-°C) 

The  value  for  empirical  coefficient  6  can  be  taken  as  0.4  (l.e.,  von 
Karmen' s  constant)  although  Sundaram  (Mortimer  1974)  recommends  0.2  to 
0.4.  The  surface  heat  flux  can  be  estimated  from 

H  *  K(T  -  T  )  (A4) 

n  e  s 

where 

2 

K  »  heat  exchange  coefficient,  W/  (m  -  °C) 

T  =  equilibrium  temperature,  °C 
T^  ■»  surface  temperature,  °C 

Procedures  for  computing  K  and  Tg  can  be  found  in  Edinger,  Brady, 
and  Geyer  (1974). 

The  physical  significance  of  is  that  it  is  a  measure  of  the 

depth  the  wind  can  distribute  a  given  surface  heat  input.  Lakes  with 
depths  greater  than  D^,  and  not  dominated  by  advection  (flowthrough) 
will  probably  stratify. 
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A.  3. 2  Examples 

Since  the  wind  mixing  depth  scale  (Equation  A3)  is  dependent  only 
on  meteorological  conditions,  computations  will  be  made  for  southern 
Minnesota  (Minneapolis)  and  central  Louisiana  (Monroe)  to  illustrate 
geographic  variations  in  .  Equation  A3  reduces  to 

5.9  x  LO9  v* 

°t - H -  <A5> 

n 

with  substitution  for  the  various  coefficients  and  constants.  Assuming 

_3 

a  drag  coefficient  of  1,3  *  10  (see  Equations  54  and  55),  the 

shear  velocity  w*  can  be  obtained  directly  from  the  wind  speed  W 
using 


1.27  x  l(f3  W 


(A6) 


where  W  ■  wind  speed,  m/sec. 

The  surface  heat  is  obtained  from  Equation  A4  using  Figure  A1  to 

obtain  K  and  approximating  T  by 
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T  +  — 
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where 


H 

s 


dew  point  temperature,  °C 

2 

gross  rate  of  solar  radiation,  W/m 


Computation  of  therefore  requires  values  for  the  wind  speed,  dew 

point  temperature,  water  surface  temperature,  and  solar  radiation  for 
the  period  of  interest.  Assuming  the  period  of  interest  for  Minnesota 
was  May  through  September  and  for  Louisiana  was  April  through  November, 
the  following  meteorological  parameters  were  obtained  from  historical 
records  and  used  as  input. 
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AVERAGE  TEMPERATURE,  Tm  *  (  Td +TS )/2  (C) 


Figure  Al.  Design  curve  for  surface  heat  exchange  coefficient  K 
(after  Edinger,  Brady,  and  Geyer  1974) 


Minnesota 


Wind  speed,  .n/sec  4.3 

Dew  point  temperature,  °C  11.3 

2 

Solar  radiation,  W/m  230 

Water  surface  temperature,  °C  15,2 


Louisiana 

2.6 

16.7 

223 

23.0 


The  computed  quantities  are: 


K  from  Figure  Al 
from  Equation  A7 
from  Equation  A4 
wA  from  Equation  A6 
from  Equation  A5 


Minnesota 

Louisiana 

29 

29 

19.2 

24.8 

117 

49.8 

0.0055 

0.0033 

8.4 

4.3 
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Since  these  Dt  values  are  less  than  the  general  rule  of  10  m,  they 
Indicate  the  importance  of  other  mixing  processes  in  determining  the 
stratification  potential. 


A. 4  Pond  Number 


A.  4 . 1  Description 

For  reregulation  pools  or  multiple  reservoirs  in  series,  a  more 
sophisticated  stratification  criterion  has  been  proposed  by  Jirka  and 
Harleman  (1979)  and  Jirka  and  Watanabe  (1980).  Although  this  criterion 
was  originally  developed  for  cooling  ponds,  it  can  be  used  on  any  system 
characterized  by  large  unsteady  inflows  of  different  density  (tempera¬ 
ture)  water,  Jirka  and  Harleman  (1979)  proposed  the  Pond  number,  Po, 
which  is  defined  by: 


Po 


f  2 

2±  Q  d3  L_ 

4  -.302  v  D 

aATgD  B  m 

m  c 


1/4 


(A8) 


where 

f^  *  interfacial  friction  coefficient  (~0.01),  dimensionless 

a  ■  coefficient  of  thermal  expansion,  per  °C 

AT  -  temperature  differential  between  inflow  water  and  ambient 
reservoir  water,  °C 

»  dilution  ratio  for  entrance  mixing  (~1.0),  dimensionless 

The  Pond  number  Includes  four  separate  factors: 

a.  The  parameter  f^/ 4  describes  the  magnitude  of  the  internal 
turbulent  shear, 

2  3  2  -1 

b.  The  densiraetric  Froude  number  Q  (aATgD  B  )  compares  the 

—  me 

destabilizing  kinetic  energy  of  inflow  with  the  stabilizing 

buoyant  energy. 

3 

c.  The  factor  considers  the  destabilizing  entrance  mixing. 

d.  The  parameter  L/D  is  an  aspect  ratio  comparing  length  to 
depth  dimensions. 
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The  larger  the  values  of  Po  ,  the  weaker  the  stratification.  If 
Po  <  0.3  ,  the  pool  is  well  stratified.  If  0.3  £  Po  £  1,0  ,  the  pool 
is  weakly  stratified  with  a  vertical  temperature  difference  of  AT  ■ 
0.45To  (1-Po)  ,  where  Tq  ■  epilimnetic  temperature,  °C.  If  Po  >  1.0  , 
the  pool  is  vertically  well  mixed. 

A. 4. 2  Examples 

Lakes  Catherine  and  Hamilton  are  two  small  hydropower  projects 
located  downstream  of  Lake  Ouachita,  a  large  CE  reservoir,  on  the 
Ouachita  River  in  south-central  Arkansas.  All  of  these  projects  have 
bottom  releases  and  are  operated  in  series  such  that  the  cold, 
hypolimnetic  releases  from  Lake  Ouachita  pass  through  Lake  Hamilton  as 
an  underflow  and  then  through  Lake  Catherine  as  an  underflow.  The  Pond 
number  criterion  is  used  to  evaluate  the  stratification  potential  of 
Lakes  Catherine  and  Hamilton  because  there  is  a  large  unnatural  tempera¬ 
ture  difference  between  the  inflowing  release  waters  from  the  upstream 
reservoir  and  natural  reservoir  surface  temperatures. 

In  this  example,  the  following  coefficients  were  used 


f±  -  0.01 

a  -  1.8  x  10"Vc 
2 

g  *  9.8  m/sec 
D  -  1.0 


reducing  Equation  A8  to 


Po  *  1.09 


,  IM 


ATD3B2  ^m 
m 


The  lakes  have  the  following  morphometric  characteristics: 


(A9) 


Characteristic,  m 

Catherine 

Hamilton 

D 

5.5 

8.0 

m 

B 

400 

970 

L 

19,500 

29,900 

A8 

and  mean  monthly  values  were  used  for  the  flows  Q  and  the  temperature 
differentials.  The  temperature  differentials  were  obtained  from  sine 
curves  fitted  to  field  measurements  (see  Figure  A2) .  The  following 
values  were  used 


Catherine 

Hamilton 

3 

May  -  Q,  m  / sec 

42.5 

42.0 

AT,  °C 

9.0 

9.8 

3 

July  -  Q,  m  /sec 

28.3 

25.0 

AT,  “C 

10.5 

10.3 

3 

October  -  Q,  m  /sec 

34.0 

27.0 

AT,  °C 

2.3 

2.5 

which,  when  substituted  into  Equation  A9,  resulted  in  the  following  Pond 
numbers : 


Catherine  Hamilton 


May 

0.44 

0.21 

July 

0.34 

0.19 

October 

0.55 

0.25 

These  results  indicate  Lake  Hamilton  should  be  well  stratified  since 
Po  <  0.3  and  Lake  Catherine  should  be  weakly  stratified  since  0,3  £  Po 
<  1.0  .  Field  measurements  verify  these  results  with  midsummer 
temperatures  in  Lake  Hamilton  varying  from  14°  C  in  the  hypolimnion  to 
30°  C  in  the  epilimnion  and  in  Lake  Catherine  varying  from  18°  C  in  the 
hypolimnion  to  28°  C  in  the  epilimnion.  In  Lake  Catherine,  the  vertical 
temperature  differential  can  be  estimated  from 

AT  =  0.45T  (1-Po) 

o 

If  Po  =  0,34  and  T  =  28°  C,  then  AT  =  8.3°  C,  which  is  similar  to 

o 

the  measured  10°  C. 
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Figure  A2.  Sine  curves  fit  to  surface,  outflow,  and  inflow 
temperatures  for  Lake  Hamilton 


A . 5  Summary 

The  methods  presented  to  evaluate  the  stratification  potential  of  a 
reservoir  range  from  a  general  rule  to  simple  computations.  Used  inde¬ 
pendently,  these  methods  can  provide  an  indication  of  the  stratification 
potential.  However,  because  of  the  many  factors  involved  in  the 
development  and  maintenance  of  stratification,  the  use  of  more  than  one 
of  the  methods  is  suggested  to  provide  a  broader  basis  for  the 
evaluation  of  the  stratification  potential  of  a  given  reservoir  (see 
Section  A.2.2) . 
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NOTATION 


A 

s 


A(z) 

a 


B 

B(z) 

C 


c 


GSWH 


g 

H 

H 

n 

H 

s 


2 

Cross-sectional  area  of  inflow,  m 
Water  surface  area,  m^ 

2 

Horizontal  area  of  reservoir  at  elevation  z,  m 

Dimensionless  coefficient 

Width  of  zone  of  conveyance,  m 

Reservoir  width  at  elevation  z,  tn 

3 

Concentration,  g/m 

Drag  coefficient,  dimensionless 

Empirical  coefficient  or  proportionality  constant,  m 
Dimensionless  calibration  coefficient 
Calibration  coefficient  for  advective  mixing 
Empirical  calibration  coefficient 
Heat  capacity  of  water,  J/(kg-°C) 

Constants  determined  from  experimental  data 
Calibration  coefficient  for  wind  mixing 
c/a  ,  m 

Water  depth,  m 
Thickness  of  top  layer,  m 
Reservoir  mean  depth,  m 
Wind  mixing  depth  scale,  m 

Dilution  ratio  for  entrance  mixing  (~I.O),  dimensionless 

Slope  of  water  surface,  m/m 

Dimensionless  density  gradient  =  10  ^/m 
Densimetric  Froude  number,  dimensionless 
Coriolis  parameter  (f  =  2W  sin  0) 

Interfacial  friction  coefficient  (~0.0i),  dimensionless 
WQRRS  stability  criterion 

2 

Acceleration  due  to  gravity,  m/sec 
Thickness,  depth,  or  height,  m 

2 

Net  heat  flux  across  the  air/water  interface,  W/m 

2 

Surface  heat  flux,  W/m 
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H 

w 

H 

z 

h 


zo 


k 

L 

M 


tn 


n 

Po 

Pr 


PE 

P? 

Q 


^tn 

qt 

Ri 


RPE 


Wave  height,  m 

2 

Heat  flux  at  elevation  z,  W/m 

Depth  (thickness)  of  upper  turbulent  layer  or  mixed 
layer,  m 

Depth  of  center  of  mass  of  the  mixed  layer,  m 
Hydraulic  plunge  depth,  m 

2 

Heat  exchange  coefficient,  W/ (tn  -°C) 

Kinetic  energy,  J 

2 

Dispersion  coefficient,  tn  /sec 

2 

Molecular  diffusivity  or  diffusion  coefficient,  m  / sec 

2 

Turbulent  diffusion  coefficient,  m  /sec 

2 

Global  vertical  diffusion  coefficient,  m  /sec 

2 

Eddy  diffusion  coefficient  at  neutral  stability,  m  /sec 
von  Kantian*  s  constant  (~0.4) 

Length  scale  or  reservoir  length,  m 
Mass,  kg 

Total  mass  of  the  reservoir,  kg 


^  4^  *  buoyancy  frequency,  sec  2 

P  az 

Coefficient,  dimensionless 
Pond  number,  dimensionless 

Prandtl  number  =  1  for  water,  dimensionless 
Potential  energy,  J 

Turbulent  fluctuations  of  pressure.  Pa 

3 

Average  reservoir  discharge,  m  /sec 

3 

Flow  rate  in  layer  i,  m  /sec 

3 

Inflow  rate,  m  /sec 

3 

Vertical  flow  rate,  m  /sec 


2  2  2 

TKE  per  unit  mass  =(u  +  u  +  u  )/2 

x  y  2  ^ 

Flux  due  to  dispersion,  kg/(m  -sec) 

2 

Flux  due  to  molecular  diffusion,  kg/(m  -sec) 

2 

Flux  due  to  turbulent  diffusion,  kg/m  -sec) 
Richardson  number,  dimensionless 
Relative  potential  energy,  J 
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r 

S 


T 


TKE 

t 

TKE 

c 

TKE 


s 

t 


s 

U 

U 


uj 

u, 


u  ,  u  ,  u 

x  y  z 


V 


u  ,  u 
y  z 

u.(z) 

uo(z) 

u 


u'w' 


u  ,  u  ,  u 
x  y  z 


u  u 
X  z 


u 


'i 

W 


Radius  of  an  inertial  circle,  m 
Water  surface  slope,  m/m 
Stability  factor,  dimensionless 
Interface  slope,  m/m 
Water  temperature,  °C 

2  2 

Turbulent  kinetic  energy  from  advection,  (kg-m  )/sec 

2  2 

Turbulent  kinetic  energy  from  convection,  (kg-m  )/sec 

2  2 

Wind  shear  turbulent  kinetic  energy,  (kg-m  )/sec 

Dew  point  temperature,  °C 

Equilibrium  temperature,  °C 

Inflow  temperature,  °C 

Surface  temperature,  °C 

Time ,  sec 

Hydraulic  residence  time,  sec 
Seiche  period,  sec 

Velocity  in  direction  of  flow,  m/sec 
Mean  horizontal  component,  m/sec 
Flow  velocity  in  layer  i,  m/sec 
Inflow  velocity,  m/sec 
Instantaneous  velocity  components,  m/s ec 

Mean  velocity  components,  m/sec 

Inflow  velocity  distribution,  m/sec 

Outflow  velocity  distribution,  m/sec 

Surface  drift  velocity,  m/sec 

2  2 

Reynolds  stress,  m  /sec 

Turbulent  fluctuations  of  the  horizontal  (x,y)  and 
vertical  (z)  velocity  components,  m/sec 

2  2 

Reynolds  stress,  m  /sec 

2  2  2  ~ 

u  +  u  +  u 
x  y  z 

Volume ,  m^ 

3 

Volume  of  layer  i,  m 
Wind  speed,  m/sec 
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2  1/2 
<u  > 


p(z,t) 


Wedderburn  number 

2  2 

Entrainment  work,  (kg-m  )/sec 
^  ■  entrainment  velocity,  m/sec 
Shear  or  friction  velocity,  m/sec 

Difference  in  depth  between  center  of  mass  of  inflow  and 
center  of  mass  of  inflow  placement,  m 

Maximum  elevation  or  depth,  m 

Vertical  coordinate,  m 

Elevation  of  water  surface,  m 

Ensemble  mean  (i.e.,  mean  over  many  trials) 

Root  mean  square  velocity,  m/sec 

Coefficient  of  thermal  expansion  of  water,  per  °C 

Coefficient,  dimensionless 

Setup  of  water  surface,  m 

Temperature  differential  between  inflow  water  and  ambient 
reservoir  water,  °C 

3 

Incremental  volume  to  be  entrained,  m 

Time  step,  sec 

Density  difference,  kg/mf' 

Density  difference  between  inflow  and  reservoir  surface 
waters,  kg/m^ 

2  3 

Rate  of  dissipation  of  TKE,  m  /sec 

2  3 

Local  dissipation  rate  of  TKE,  m  /sec 
Wavelength,  m 

3 

Density,  kg/m 

2 

Turbulent  fluctuations  of  density,  kg/m 

3 

Density  of  air,  1.177  kg/m 

3 

Inflow  density,  kg/m 

3 

Maximum  water  density,  kg/m 

3 

Mean  density,  kg/m 

3 

Water  density,  kg/m 

3 

Reservoir  density  at  elevation  z  and  time  t,  kg/m 
Velocity  scale,  m/sec 
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2 

Shear  stress,  kg/(m-sec  ) 

2 

Shear  stress  at  air/water  interface,  kg/(m-sec  ) 
Latitude 

Concentration  gradient,  mg/ i/m 
Velocity  gradient,  (m/sec)/m 
Local  density  gradient 

Angular  velocity  of  the  earth's  rotation,  7.29  x  10 
rad/sec 
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